Skip to content
GitLab
Commit 5b6d3757 authored by Luca Mastrolorenzo's avatar Luca Mastrolorenzo
Browse files

added GammaShowerInfo.hh

parent 6b61388d
Hide whitespace changes
Inline Side-by-side
interface/GammaShowerInfo.hh 0 → 100644
View file @ 5b6d3757
#ifndef GammaShowerInfo_HH
#define GammaShowerInfo_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"
//#include "coeffiecientA.hh"
double EnergyCutxLayerPhoton = 0.;
class GammaShowerInfo : public AnalysisBase {
public:
enum {
kCutBins=50
};
GammaShowerInfo(const std::string &sName="GammaShowerInfo") :
AnalysisBase(sName), fTwoPi(2.0*acos(-1.0)) {
h_PhotonE = new TH1F("h_PhotonE","Photon Energy Distribution",100, 0,100);
h_PhotonPt = new TH1F("h_PhotonPt","Photon Pt Distribution",100, 0,100);
h_PhotonEta = new TH1F("h_PhotonEta","Photon #eta Distribution",100, -4, 4);
h_PhotonPhi = new TH1F("h_PhotonEPhi","Photon #phi Distribution",100, -4, 4);
//h_Matrix = new TH2F("h_Matrix","d_e_i * d_e_j", 40, 0, 40, 40, 0, 40);
h_Matrix = new TH2F("h_Matrix","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 = new TH1F("h_Vector","vector", 53, 0, 53);
h_N = new TH1F("h_N","N", 53, 0, 53);
h_mipInWafer = new TH1F("h_mipInWafer","h_mipInWafer", 250, 0, 1000);
h_mipInLayer = new TH1F("h_mipInLayer","h_mipInLayer", 53, 0, 53);
h_mipInEndcap = new TH1F("h_mipInEndcap","h_mipInEndcap", 250, 0, 50000);
h_mipVsLayer = new TH2F("h_mipVsLayer","MIP vs Layer", 53, 0, 53, 500, 0, 500);
h_EtInLayer = new TH1F("h_EtInLayer","h_EtInLayer", 53, 0, 53);
h_EcalibVsEta = new TH2F("h_EcalibVsEta","h_EcalibVsEta", 500, 0, 500, 100, 0, 3.1);
h_MaxEposition = new TH1F("h_MaxEposition","h_MaxEposition", 53, 0, 53);
h_MaxEpositionVsEnergy = new TH2F("h_MaxEpositionVsEnergy","h_MaxEpositionVsEnergy", 53, 0, 53, 250, 0, 5000);
h_PiReso_step0 = new TH1F("h_PiReso_step0","h_PiReso_step0", 100, 0, 2);
h_PiReso_step0_Et = new TH1F("h_PiReso_step0_Et","h_PiReso_step0_Et", 100, 0, 2);
h_PiResoCalib_step0 = new TH1F("h_PiResoCalib_step0","h_PiResoCalib_step0", 100, 0, 2);
h_EfirstLayer_vs_R = new TH2F("h_EfirstLayer_vs_R","h_EfirstLayer_vs_R", 100, -2, 2, 500, 0, 500);
for(unsigned l(0);l<=Geometry::kNumberOfLayers;l++) {
h_ELayer_vs_R[l]=new TH2F((sName+sLayerLabel[l]+"h_ELayer_vs_R").c_str(),
(sLayerTitle[l]+";MIP vs Resolution").c_str(),
100, -2, 2, 500, 0, 500);
}
for(unsigned l(0);l<=Geometry::kNumberOfLayers;l++) {
h_HitFracXLayer[l]=new TH1F((sName+sLayerLabel[l]+"h_HitFracXLayer").c_str(),
(sLayerTitle[l]+";Fraction of the recovered hits").c_str(),
100,0.0,1.2);
}
//double calib( double dep[] );
}
virtual ~GammaShowerInfo() {
}
virtual bool event(Event &event) {
bool debug = false;
std::cout << fName << " Event " << fEventNumber << std::endl;
for(unsigned e(0);e<Geometry::kNumberOfEndcaps;e++) {
double t_deposits[Geometry::kNumberOfLayers + Geometry::kNumberOfBHLayers+1];
double t_weightDep[Geometry::kNumberOfLayers + Geometry::kNumberOfBHLayers+1];
std::vector<Particle> &vPhoton(event.photons(e));
bool has_lzero_2mips = false;
//if(!( fabs(vPi[0].FourP().Eta())>2.0 && fabs(vPi[0].FourP().Eta())<2.1 ) ) continue;
if(debug)std::cout << "e = " << e << " photon size " << vPhoton.size() <<std::endl;
double Et_inPiCone=0.;
for(unsigned i(0);i<vPhoton.size();i++) {
h_PhotonE->Fill(vPhoton[i].FourP().E());
h_PhotonPt->Fill(vPhoton[i].FourP().Pt());
h_PhotonEta->Fill(vPhoton[i].FourP().Eta());
h_PhotonPhi->Fill(vPhoton[i].FourP().Phi());
if(debug)std::cout << "photon " << i << " (E,pt,eta,phi) = " << vPhoton[i].FourP().E() << ", " << vPhoton[i].FourP().Pt() << ", " << vPhoton[i].FourP().Eta() << ", " << vPhoton[i].FourP().Phi() << std::endl;
}
double EmEt=0.;
double HadEt=0.;
double Et_TOT_endcapInCone = 0.;
double MIP_TOT_endcap = 0.;
double Et_TOT_endcap = 0.;
double E_TOT_endcap = 0.;
bool IsOver[Geometry::kNumberOfLayers];
bool useful(false);
//int useful = 0;
int goodLayer = 0;
for(unsigned l(0); l<Geometry::kNumberOfLayers; l++) {
std::vector<Particle> vPhotonProjTest;
for(unsigned i(0);i<vPhoton.size();i++){
vPhotonProjTest.push_back(vPhoton[i].project(Geometry::layerZ(e,l)));
}
for(unsigned i(0);i<vPhotonProjTest.size();i++){
if( vPhotonProjTest.at(i).position().rho() > (31.9 + 7) && vPhotonProjTest.at(i).position().rho() < (150 - 7) ) goodLayer = goodLayer + 1 ;
}
}
if(goodLayer<40) continue;
double tmp_ELayer[Geometry::kNumberOfLayers];
for(unsigned l(0);l<Geometry::kNumberOfLayers;l++) {
//std::cout << "LAYER " << l << std::endl;
// Find projection of pion position into layer, and fill a vector at each layer
std::vector<Particle> vPhotonProj;
//loop over pion and assign for each layer its corresponding projection
for(unsigned i(0);i<vPhoton.size();i++){
vPhotonProj.push_back(vPhoton[i].project(Geometry::layerZ(e,l)));
}
double tMIP_TOT_layer[100];
double tMIP_TOT_layerInCone[100];
double tEt_TOT_layerInCone[100];
//look to the pion projection in each layer and see if they satisfy the requirements
if(debug)std::cout << "projection size = "<<vPhotonProj.size() << std::endl;
//if( ! ( vPhotonProj[0].position().rho() > 31.9 + 5 && vPhotonProj[39].position().rho() < 150 - 5) ) continue;
for(unsigned i(0);i<vPhotonProj.size();i++){
if(debug)std::cout << " ----PiProj Pt "<< vPhotonProj[i].FourP().Eta() << " rho = " << vPhotonProj[i].position().rho() << std::endl;
if( !(vPhotonProj.size()>0 && vPhotonProj[i].FourP().Pt()>0 && vPhotonProj[i].position().rho()<130.0 && vPhotonProj[i].position().rho()>40.0 ) ) continue;
useful=true;
double MIP_TOT_layer = 0.;
double MIP_TOT_layerInCone = 0.;
double Et_TOT_layer = 0.;
double Et_TOT_layerInCone = 0.;
double E_TOT_layer = 0.;
for(unsigned w(0);w<Geometry::numberOfWafers(l);w++) {
if(Geometry::validWafer(l,w)) {
double MIP_TOT_wafer = 0.;
double Et_TOT_wafer = 0.;
double E_TOT_wafer = 0.;
//loop over the cells of that wafer (in that layer)
for(unsigned c(0);c<Geometry::numberOfCells(l,w);c++) {
//retrieve the point of the cell and the vector of the simHit that have hit the cell
Point p(Geometry::waferCellPoint(e,l,w,c));
std::vector<SimHit> &vSimHit(event.simHits(e,l,w,c));
double sumCell(0.0);
double sumCellEt(0.0);
double sumCellE(0.0);
//loop over the simHit of the cell and sum all their energy ==> obtaining the total energy deposited in the cell (in MIP)
for(unsigned s(0);s<vSimHit.size();s++){
sumCell+=vSimHit[s].mips();
//sumCell+=vSimHit[s].deposit();
sumCellEt+=vSimHit[s].transverseEnergy();
sumCellE+=vSimHit[s].energy();
}
if(sumCell>0){if(debug)std::cout << "energy in the cell = "<< l << " , " << w << " , " << c << ", " << sumCell << std::endl;
MIP_TOT_wafer += sumCell;
Et_TOT_wafer += sumCellEt;
E_TOT_wafer += sumCellE;
}
}
if(MIP_TOT_wafer>0){if(debug)std::cout << " energy in the wafer = " << l << " , " << w << " , " << MIP_TOT_wafer << std::endl;
MIP_TOT_layer += MIP_TOT_wafer;
Et_TOT_layer += Et_TOT_wafer;
E_TOT_layer += E_TOT_wafer;
h_mipInWafer->Fill( MIP_TOT_wafer );
}
}
}
// if((MIP_TOT_layer>0 && MIP_TOT_layer<2 && l==0) || ( MIP_TOT_layer>0 && (l!=0 || l!=39) ) || (MIP_TOT_layer>0 && MIP_TOT_layer<2 && l==39) ){
tmp_ELayer[l]=MIP_TOT_layer;
//if(l==0 && MIP_TOT_layer>4)
// has_lzero_2mips = true;
//if(has_lzero_2mips)
// continue;
if(MIP_TOT_layer>EnergyCutxLayerPhoton) IsOver[l]=true;
else IsOver[l]=false;
if(MIP_TOT_layer>0){
if(debug)std::cout << " energy in the layer = " << l << " , " << MIP_TOT_layer << std::endl;
h_mipInLayer->Fill(l, MIP_TOT_layer );
h_mipVsLayer->Fill(l, MIP_TOT_layer );
MIP_TOT_endcap += MIP_TOT_layer;
Et_TOT_endcap += Et_TOT_layer;
E_TOT_endcap += E_TOT_layer;
//tMIP_TOT_layer[i] = MIP_TOT_layer;
//tMIP_TOT_layerInCone[i] = MIP_TOT_layerInCone;
if( MIP_TOT_layer >0 ){
//std::cout << "Layer "<< l << " total energy = " << tMIP_TOT_layer[i] << " total energy in the cone = " << tMIP_TOT_layerInCone[i] << " fraction of collected energy = " << tMIP_TOT_layerInCone[i] / tMIP_TOT_layer[i] << std::endl;
//h_HitFracXLayer[l]->Fill(tMIP_TOT_layerInCone[i] / tMIP_TOT_layer[i]);
//h_HitFracXLayer[Geometry::kNumberOfLayers]->Fill(tMIP_TOT_layerInCone[i] / tMIP_TOT_layer[i]);
//h_mipInLayer->Fill(l, tMIP_TOT_layer[i]);
// h_EtInLayer->Fill(l, tEt_TOT_layerInCone[i]);
t_deposits[l] = MIP_TOT_layer;
if(debug)std::cout << " ------------ deposit = " << t_deposits[l]<< std::endl;
//it works because there is one pion per endcap!!
t_weightDep[l] = MIP_TOT_layer * vPhotonProj[i].FourP().E();
}else{
t_deposits[l] = 0;
t_weightDep[l] = 0;
}
}
}
}//end loop layer
t_deposits[52] = 1.;
t_weightDep[52] = 1. * vPhoton[0].FourP().E();
std::cout << "check = " << t_deposits[52] << " , " << t_weightDep[52] << std::endl;
for(unsigned lbh(0); lbh<Geometry::kNumberOfBHLayers; lbh++) {
double MIP_TOT_layerBH(0.0);
double Et_TOT_layerBH(0.0);
double E_TOT_layerBH(0.0);
//loop to retrieve the BH
for(unsigned int h=15; h<99 ; h++ ){
for(unsigned int f=0; f<359 ; f++){
std::vector<SimHbh> &vSimHitBH(event.simHbhs(e,lbh,h,f));
for(unsigned sbh(0);sbh<vSimHitBH.size();sbh++){
MIP_TOT_layerBH+=vSimHitBH[sbh].mips();
//sumCell+=vSimHit[s].deposit();
Et_TOT_layerBH+=vSimHitBH[sbh].transverseEnergy();
E_TOT_layerBH+=vSimHitBH[sbh].energy();
}
//if(MIP_TOT_layerBH>0)std::cout << " BH = " << MIP_TOT_layerBH << std::endl;
}
}
h_mipInLayer->Fill(Geometry::kNumberOfLayers+lbh, MIP_TOT_layerBH );
h_mipVsLayer->Fill(Geometry::kNumberOfLayers+lbh, MIP_TOT_layerBH );
if( MIP_TOT_layerBH >0 ){
//std::cout << "Layer "<< l << " total energy = " << tMIP_TOT_layer[i] << " total energy in the cone = " << tMIP_TOT_layerInCone[i] << " fraction of collected energy = " << tMIP_TOT_layerInCone[i] / tMIP_TOT_layer[i] << std::endl;
//h_HitFracXLayer[l]->Fill(tMIP_TOT_layerInCone[i] / tMIP_TOT_layer[i]);
//h_HitFracXLayer[Geometry::kNumberOfLayers]->Fill(tMIP_TOT_layerInCone[i] / tMIP_TOT_layer[i]);
//h_mipInLayer->Fill(l, tMIP_TOT_layer[i]);
// h_EtInLayer->Fill(l, tEt_TOT_layerInCone[i]);
t_deposits[Geometry::kNumberOfLayers+lbh] = MIP_TOT_layerBH;
if(debug)std::cout << " ------------ deposit = " << t_deposits[Geometry::kNumberOfLayers + lbh]<< std::endl;
//it works because there is one pion per endcap!!
t_weightDep[Geometry::kNumberOfLayers+lbh] = MIP_TOT_layerBH * vPhoton[0].FourP().E();
}else{
t_deposits[Geometry::kNumberOfLayers+lbh] = 0;
t_weightDep[Geometry::kNumberOfLayers+lbh] = 0;
}
}
int nLayerGood=0;
for(int hi=0; hi<Geometry::kNumberOfLayers; hi++){
if(hi<10) continue;
if(IsOver[hi]) nLayerGood++;
}
// if(MIP_TOT_endcap>0 && nLayerGood>3){
if(MIP_TOT_endcap>0){
if(debug)std::cout << " energy in the ENDCAP = " << e << " , " << MIP_TOT_endcap << std::endl;
h_mipInEndcap->Fill( MIP_TOT_endcap );
h_PiReso_step0->Fill( E_TOT_endcap / vPhoton[0].FourP().E() );
h_PiReso_step0_Et->Fill( Et_TOT_endcap / vPhoton[0].FourP().Et() );
double E_PiCalib = calib(t_deposits);
h_PiResoCalib_step0->Fill( E_PiCalib / vPhoton[0].FourP().E() );
h_EfirstLayer_vs_R->Fill( E_PiCalib / vPhoton[0].FourP().E() ,tmp_ELayer[0] );
for(int kk=0; kk<40; kk++){
h_ELayer_vs_R[kk]->Fill( E_PiCalib / vPhoton[0].FourP().E() ,tmp_ELayer[kk] );
}
h_EcalibVsEta->Fill( E_PiCalib, fabs( vPhoton[0].FourP().Eta() ) );
}
for(int i=0; i<Geometry::kNumberOfLayers+Geometry::kNumberOfBHLayers+1; i++){
h_N->Fill(i,1);
if(t_weightDep[i]>0){
h_Vector->Fill(i,t_weightDep[i]);
if(debug)std::cout << " =============== fill vector = " << t_weightDep[i] << std::endl;
}
for(int j=0; j<Geometry::kNumberOfLayers+Geometry::kNumberOfBHLayers+1; j++){
h_N2->Fill(i, j, 1);
//std::cout <<h_N2->GetBinContent(i,j) << std::endl;
if(t_deposits[i]>0 && t_deposits[j]>0){
//std::cout << "i,j " << i << ", "<< j<< " dep_i * dep_j = " << t_deposits[i] * t_deposits[j] <<std::endl;
//h_Matrix->Fill( i, j, t_deposits[i] * t_deposits[j] );
h_Matrix->Fill( i, j, t_deposits[i] * t_deposits[j] );
//std::cout << " AFTER i,j " << i << ", "<< j<< " dep_i * dep_j = " << t_deposits[i] * t_deposits[j] <<std::endl;
}
}
}
int MaxEposition = distance( t_deposits, std::max_element( t_deposits, t_deposits+53 ) );
std::cout << "max element" << max_element( t_deposits, t_deposits+53 ) << " index" << MaxEposition<< std::endl;
double MaxE = t_deposits[MaxEposition];
h_MaxEposition->Fill(MaxEposition);
h_MaxEpositionVsEnergy->Fill(MaxEposition, MaxE);
}
return true;
}
private:
const double fTwoPi;
TH2F *h_Matrix;
TH1F *h_HitFracXLayer[Geometry::kNumberOfLayers+1];
TH2F *h_ELayer_vs_R[Geometry::kNumberOfLayers+1];
//TH2F *h_Matrix;
TH2F *h_N2;
TH1F *h_Vector;
TH1F *h_N;
TH1F *h_mipInWafer;
TH1F *h_mipInLayer;
TH1F *h_mipInEndcap;
TH1F *h_EtInLayer;
TH2F *h_mipVsLayer;
TH1F *h_MaxEposition;
TH2F *h_MaxEpositionVsEnergy;
TH1F *h_PiReso_step0;
TH1F *h_PiReso_step0_Et;
TH1F *h_PiResoCalib_step0;
TH2F *h_EfirstLayer_vs_R;
TH2F *h_EcalibVsEta;
TH1F *h_PhotonE;
TH1F *h_PhotonPt;
TH1F *h_PhotonEta;
TH1F *h_PhotonPhi;
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;
}
};
#endif
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment