Skip to content

GitLab

Commit e91c0f5f authored by Luca Mastrolorenzo's avatar Luca Mastrolorenzo
Browse files
Options

Added new class to analyze cluster-3D CoefficientDetClu3D.hh

parent 6a80a1fa
Hide whitespace changes
Inline Side-by-side
Showing with 388 additions and 3 deletions
interface/CoefficientDetClu3D.hh 0 → 100644
View file @ e91c0f5f
#ifndef CoefficientDetClu3D_HH
#define CoefficientDetClu3D_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 "coefficientCluster.hh"
//#include "Calibration.hh"
double Eta_min = 1.7;
double Eta_max = 2.7;
TString ParticleType = "photon";
class CoefficientDetClu3D : public AnalysisBase {
public:
enum {
kCutBins=50
};
CoefficientDetClu3D(const std::string &sRoot="", const std::string &sOut="")
: AnalysisBase("CoefficientDetClu3D",sRoot,sOut) {
//AnalysisBase(sName), fTwoPi(2.0*acos(-1.0)) {
h_SigE = new TH1F("h_SigE","Signal Energy Distribution",200, 0,500);
h_SigPt = new TH1F("h_SigPt","Signal Pt Distribution",100, 0,200);
h_SigEta = new TH1F("h_SigEta","Signal #eta Distribution",100, -4, 4);
h_SigPhi = new TH1F("h_SigEPhi","Signal #phi Distribution",100, -4, 4);
h_FracMipTMaxC3d = new TH1F("h_FracMipTMaxC3d","h_FracMipTMaxC3d",110, 0, 1.1);
h_layer_MipTweight = new TH1F("h_layer_MipTweight","h_layer_MipTweight", 52,0,52);
h_Matrix = new TH2F("h_Matrix","d_e_i * d_e_j", 53, 0, 53, 53, 0, 53);
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 = new TH1F("h_Vector","vector", 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_CluResponse = new TH1F("h_CluResponse","h_CluResponse",100,0,2);
h_CluResponse2D = new TH2F("h_CluResponse2D","h_CluResponse2D",600,0,300,600,0,300);
// h_CluEtaResponse = new TH1F("h_CluEtaResponse", "h_CluEtaResponse", 200, -0.5, 0.5);
// h_CluPhiResponse = new TH1F("h_CluPhiResponse", "h_CluPhiResponse", 200, -0.5, 0.5);
// h_CluEtaResponseVsEtaTrue = new TH2F("h_CluEtaResponseVsEtaTrue","h_CluEtaResponseVsEtaTrue",200,1.4,3.2, 200, -0.5, 0.5);
// h_CluPhiResponseVsPhiTrue = new TH2F("h_CluPhiResponseVsPhiTrue","h_CluPhiResponseVsPhiTrue",200,-3.5,3.5, 200, -0.5, 0.5);
h_CluEtaResponseVsEtaC2d = new TH2F("h_CluEtaResponseVsEtaC2d","h_CluEtaResponseVsEtaC2d",200,1.4,3.2, 200, -0.5, 0.5);
h_CluPhiResponseVsPhiC2d = new TH2F("h_CluPhiResponseVsPhiC2d","h_CluPhiResponseVsPhiC2d",200,-3.5,3.5, 200, -0.5, 0.5);
h_CluEtaResponseVsEtTrue = new TH2F("h_CluEtaResponseVsEtTrue","h_CluEtaResponseVsEtTrue",150,0,150, 200, -0.5, 0.5);
h_CluPhiResponseVsEtTrue = new TH2F("h_CluPhiResponseVsEtTrue","h_CluPhiResponseVsEtTrue",150,0,150, 200, -0.5, 0.5);
h_CluEtaC2dVsEtaTrue = new TH2F("h_CluEtaC2dVsEtaTrue","h_CluEtaC2dVsEtaTrue", 200, 1.4, 3.2, 200, 1.4, 3.2);
h_CluPhiC2dVsPhiTrue = new TH2F("h_CluPhiC2dVsPhiTrue","h_CluPhiC2dVsPhiTrue", 200, -3.5, 3.5, 200, -3.5, 3.5);
for(unsigned e(0); e<=Geometry::kNumberOfEndcaps; e++){
std::cout << ("endcap name "+sEndcapLabel[e]).c_str() << std::endl;
for(unsigned l(0); l<=Geometry::kNumberOfLayers; l++){
h_CluEtaResponse[e][l] = new TH1F(("h_CluEtaResponse_"+sEndcapLabel[e]+"_"+sLayerLabel[l]).c_str(), ("h_CluEtaResponse_"+sEndcapTitle[e]+"_"+sLayerTitle[l]).c_str(), 200, -0.5, 0.5);
h_CluPhiResponse[e][l] = new TH1F(("h_CluPhiResponse_"+sEndcapLabel[e]+"_"+sLayerLabel[l]).c_str(), ("h_CluPhiResponse_"+sEndcapTitle[e]+"_"+sLayerTitle[l]).c_str(), 200, -0.5, 0.5);
h_CluEtaResponseVsEtaTrue[e][l] = new TH2F(("h_CluEtaResponseVsEtaTrue_"+sEndcapLabel[e]+"_"+sLayerLabel[l]).c_str(),("h_CluEtaResponseVsEtaTrue_"+sEndcapTitle[e]+"_"+sLayerTitle[l]).c_str(),200,1.4,3.2, 200, -0.5, 0.5);
h_CluPhiResponseVsPhiTrue[e][l] = new TH2F(("h_CluPhiResponseVsPhiTrue_"+sEndcapLabel[e]+"_"+sLayerLabel[l]).c_str(),("h_CluPhiResponseVsPhiTrue_"+sEndcapTitle[e]+"_"+sLayerTitle[l]).c_str(),200,-3.5,3.5, 200, -0.5, 0.5);
}
}
h_CluXResponse = new TH1F("h_CluXResponse", "h_CluXResponse", 200, -10, 10);
h_CluYResponse = new TH1F("h_CluYResponse", "h_CluYResponse", 200, -10, 10);
h_SimMipTmaxC2d_SimMipTLayer = new TH1F("h_SimMipTmaxC2d_SimMipTLayer","h_SimMipTmaxC2d_SimMipTLayer",100,0,2);
h_SimMipTallC2d_SimMipTLayer = new TH1F("h_SimMipTallC2d_SimMipTLayer","h_SimMipTallC2d_SimMipTLayer",100,0,2);
h_MipTmaxC2d_SimMipTmaxC2d = new TH1F("h_MipTmaxC2d_SimMipTmaxC2d","h_MipTmaxC2d_SimMipTmaxC2d",100,0,2);
h_dAB = new TH2F("h_dAB","h_dAB",600, 300,360,100,0,25);
h_dAB_vs_SigEt = new TH2F("h_dAB_vs_SigEt","h_dAB_vs_SigEt",500,0,500,100,0,25);
h_RatoMipTmaxC2d_MipTsecmaxC2d = new TH1F("h_RatoMipTmaxC2d_MipTsecmaxC2d","h_RatoMipTmaxC2d_MipTsecmaxC2d",20,0,10);
h_dmaxC2dSig_vs_dsecmaxC2dSig = new TH2F("h_dmaxC2dSig_vs_dsecmaxC2dSig","h_dmaxC2dSig_vs_dsecmaxC2dSig",100,0,50,100,0,50);
h_C2dMipTdensity_vs_C2dMipT = new TH2F("h_C2dMipTdensity_vs_C2dMipT","h_C2dMipTdensity_vs_C2dMipT",100,0,200,100,0,200);
h_NTrgHitPerC2d = new TH1F("h_NTrgHitPerC2d", "h_NTrgHitPerC2d", 50, 0, 50);
h_NTrgHitPerLayer = new TH1F("h_NTrgHitPerLayer", "h_NTrgHitPerLayer", 75, 0, 75);
h_NC2dPerLayer = new TH1F("h_NC2dPerLayer", "h_NC2dPerLayer", 15, 0, 15);
h_PhiTrueFromXY = new TH1F("h_PhiTrueFromXY","h_PhiTrueFromXY",100,-3.5,7);
h_PhiTrueFrom4P = new TH1F("h_PhiTrueFrom4P","h_PhiTrueFrom4P",100,-3.5,7);
h_PtReso = new TH1F("h_PtReso","h_PtReso",100, 0, 2);
//for(unsigned l(0);l<=Geometry::kNumberOfLayers;l++) {
// h_ELayer_vs_R[l]=new TH2F((fName+sLayerLabel[l]+"h_ELayer_vs_R").c_str(),
// (sLayerTitle[l]+";MIP vs Resolution").c_str(),
// 100, -2, 2, 500, 0, 500);
//}
double deltaEta(double eta1, double eta2);
double deltaPhi(double phi1, double phi2);
double dist2D(double x1, double y1, double x2, double y2);
}
virtual ~CoefficientDetClu3D() {
}
virtual bool event(Event &event) {
if(fEventNumber <= 2){
for(int i=0; i<40; i++){
std::cout << "layer weight = " << Calibration::layerWeight(i) << std::endl;
}
}
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];
double t_deposits_T[Geometry::kNumberOfLayers + Geometry::kNumberOfBHLayers + 1];
double t_weightDep_T[Geometry::kNumberOfLayers + Geometry::kNumberOfBHLayers + 1];
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);
}
if(debug)std::cout << "e = " << e << " signal size " << vSig.size() <<std::endl;
for(unsigned i(0);i<vSig.size();i++) {
h_SigE->Fill(vSig[i].FourP().E());
h_SigPt->Fill(vSig[i].FourP().Pt());
h_SigEta->Fill(vSig[i].FourP().Eta());
h_SigPhi->Fill(vSig[i].FourP().Phi());
if(debug)std::cout << "signal " << i << " (E,pt,eta,phi) = " << vSig[i].FourP().E() << ", " << vSig[i].FourP().Pt() << ", " << vSig[i].FourP().Eta() << ", " << vSig[i].FourP().Phi() << std::endl;
}
//select event with signal particle in the geometrical and kinematical acceptance
bool useful(false);
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 ;
}
}
if(goodLayer<40) continue;
/** LOOP 3d clusters **/
int TotC2dInC3d=0;
std::vector<TrgC3d> &vTrgC3d( event.trgC3ds(e) );
std::sort( vTrgC3d.begin(), vTrgC3d.end(), eneTranOrderC3d );
//hNC3d[e]->Fill( vTrgC3d.size() );
std::cout << "------------- Number of 3D clusters: " << vTrgC3d.size() << std::endl;
double totMipT=0.;
if(vTrgC3d.size()<=0) continue;
for( std::vector<TrgC3d>::iterator clu3d = vTrgC3d.begin(); clu3d != vTrgC3d.end(); clu3d++ ) {
std::cout << "C3d mipT = " << clu3d->transverseMips() << std::endl;
totMipT += clu3d->transverseMips();
//if( clu3d->containsSigEne() ){
//
//}
/* loop over 2d clusters in 3d cluster*/
std::vector<TrgC2d*> &vTrgC2d( clu3d->trgC2ds() );
TotC2dInC3d += vTrgC2d.size();
unsigned nC2dSigEne = 0;
if( vTrgC2d.size() > 3 ){
for( int iC2d = 0; iC2d < vTrgC2d.size(); iC2d++ ) {
// std::cout << "------- C2d energy = " << vTrgC2d[iC2d]->transverseMips() << " layer " << vTrgC2d[iC2d]->layer() << std::endl;
h_layer_MipTweight->Fill(vTrgC2d[iC2d]->layer(), vTrgC2d[iC2d]->transverseMips() );
if( vTrgC2d[iC2d]->containsSigEne() )
nC2dSigEne++;
}// C2d
}
} // end loop over C3d
//select the best 3D-cluster per given endcap and assign the correspondent vector of 2D-cluster
std::vector<TrgC2d*> &vTrgC2d_fromMaxC3d( vTrgC3d[0].trgC2ds() );
std::vector<TrgC2d*> &vTrgC2d_fromSecMaxC3d( vTrgC3d[1].trgC2ds() );
if(vTrgC2d_fromMaxC3d.size()>3){
for( int iC2d = 0; iC2d < vTrgC2d_fromMaxC3d.size(); iC2d++ ) {
std::cout << "------- C2d energy = " << vTrgC2d_fromMaxC3d[iC2d]->transverseMips() << " layer " << vTrgC2d_fromMaxC3d[iC2d]->layer() << std::endl;
h_layer_MipTweight->Fill(vTrgC2d_fromMaxC3d[iC2d]->layer(), vTrgC2d_fromMaxC3d[iC2d]->transverseMips() );
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();
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] )) std::cout << "IS NAN " << std::endl;
}
}
}// C2d
}
// if(vTrgC2d_fromSecMaxC3d.size()>3){
//
// for( int iC2d = 0; iC2d < vTrgC2d_fromSecMaxC3d.size(); iC2d++ ) {
// std::cout << "------- C2d energy = " << vTrgC2d_fromSecMaxC3d[iC2d]->transverseMips() << " layer " << vTrgC2d_fromSecMaxC3d[iC2d]->layer() << std::endl;
// h_layer_MipTweight->Fill(vTrgC2d_fromSecMaxC3d[iC2d]->layer(), vTrgC2d_fromSecMaxC3d[iC2d]->transverseMips() );
// if(vTrgC2d_fromSecMaxC3d[iC2d]->transverseMips()>0 ){
// t_deposits_T[vTrgC2d_fromSecMaxC3d[iC2d]->layer()] += vTrgC2d_fromSecMaxC3d[iC2d]->transverseMips();
// t_weightDep_T[vTrgC2d_fromSecMaxC3d[iC2d]->layer()] += vTrgC2d_fromSecMaxC3d[iC2d]->transverseMips() * vSig[0].FourP().Et();
// std::cout << "--------->>>>>>>>>> " << t_deposits_T[vTrgC2d_fromSecMaxC3d[iC2d]->layer()] << std::endl;
// }else{
// t_deposits_T[vTrgC2d_fromSecMaxC3d[iC2d]->layer()] += 0;
// t_weightDep_T[vTrgC2d_fromSecMaxC3d[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] )) std::cout << "IS NAN " << std::endl;
// }
// }
// }// C2d
//
// }
t_deposits_T[52] = 1.;
t_weightDep_T[52] = 1. * vSig[0].FourP().Et();
if(true) std::cout << "check dep in l=52= " << t_deposits_T[52] << " , " << t_weightDep_T[52] << std::endl;
for(int hh=0; hh<52; hh++){
if( isnan( t_deposits_T[hh] )){
t_deposits_T[hh] = 0.;
std::cout << "IS NAN " << std::endl;
}
std::cout << "cross check deposit vec = "<< hh << ", " << t_deposits_T[hh] << std::endl;
}
double E_CluCalibMIPt_26trgLayers = calibCluster(t_deposits_T);
h_CluResponse2D->Fill(vSig[0].FourP().Et(), E_CluCalibMIPt_26trgLayers);
if(vTrgC3d.size()>0){
double fracMaxC3d = vTrgC3d[0].transverseMips() / totMipT;
std:: cout << "fraction of the energy in the main cluster: " << fracMaxC3d << std::endl;
h_FracMipTMaxC3d->Fill(fracMaxC3d);
std::cout << "energy particle mc << " << vSig[0].FourP().Et() << " energy after " << E_CluCalibMIPt_26trgLayers<< std::endl;
h_PtReso->Fill(E_CluCalibMIPt_26trgLayers/vSig[0].FourP().Et());
}else{
std:: cout << "fraction of the energy in the main cluster: " << 0. << std::endl;
h_FracMipTMaxC3d->Fill(0.);
std::cout << "energy particle mc << " << 0<< " energy after " << 0 << 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(true)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;
TH2F *h_Matrix_T;
TH2F *h_N2;
TH1F *h_Vector;
TH1F *h_Vector_T;
TH1F *h_N;
TH1F *h_SigE;
TH1F *h_SigPt;
TH1F *h_SigEta;
TH1F *h_SigPhi;
TH1F *h_CluResponse;
TH2F *h_CluResponse2D;
// TH1F *h_CluEtaResponse;
// TH1F *h_CluPhiResponse;
// TH2F *h_CluEtaResponseVsEtaTrue;
// TH2F *h_CluPhiResponseVsPhiTrue;
TH1F *h_CluEtaResponse[Geometry::kNumberOfEndcaps+1][Geometry::kNumberOfLayers+1];
TH1F *h_CluPhiResponse[Geometry::kNumberOfEndcaps+1][Geometry::kNumberOfLayers+1];
TH2F *h_CluEtaResponseVsEtaTrue[Geometry::kNumberOfEndcaps+1][Geometry::kNumberOfLayers+1];
TH2F *h_CluPhiResponseVsPhiTrue[Geometry::kNumberOfEndcaps+1][Geometry::kNumberOfLayers+1];
TH2F *h_CluEtaResponseVsEtaC2d;
TH2F *h_CluPhiResponseVsPhiC2d;
TH2F *h_CluEtaResponseVsEtTrue;
TH2F *h_CluPhiResponseVsEtTrue;
TH2F *h_CluEtaC2dVsEtaTrue;
TH2F *h_CluPhiC2dVsPhiTrue;
TH1F *h_CluXResponse;
TH1F *h_CluYResponse;
TH1F *h_SimMipTmaxC2d_SimMipTLayer;
TH1F *h_SimMipTallC2d_SimMipTLayer;
TH1F *h_MipTmaxC2d_SimMipTmaxC2d;
TH2F *h_dAB;
TH2F *h_dAB_vs_SigEt;
TH1F *h_RatoMipTmaxC2d_MipTsecmaxC2d;
TH2F *h_dmaxC2dSig_vs_dsecmaxC2dSig;
TH2F *h_C2dMipTdensity_vs_C2dMipT;
TH1F *h_NTrgHitPerC2d;
TH1F *h_NTrgHitPerLayer;
TH1F *h_NC2dPerLayer;
TH1F *h_PhiTrueFromXY;
TH1F *h_PhiTrueFrom4P;
TH1F *h_FracMipTMaxC3d;
TH1F *h_layer_MipTweight;
TH1F *h_PtReso;
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
interface/Event.hh
View file @ e91c0f5f
......@@ -53,6 +53,11 @@ bool eneTranOrderC2d( TrgC2d A, TrgC2d B)
return A.transverseMips() > B.transverseMips();
}
bool eneTranOrderC3d( TrgC3d A, TrgC3d B)
{
return A.transverseMips() > B.transverseMips();
}
class Event {
public:
......
interface/coefficientCluster.hh
View file @ e91c0f5f
......@@ -3,16 +3,20 @@ double calibCluster(double dep[]){
double E=0.;
//coefficient from pions - most energetic 2Dcluseter/layer - 1.6<|eta|<2.8 - Seed Threshold = 10 mips, TrgHit Threshold = 5 mips
double MaxCluCoeff_26l[52] = { 0, 0.0877179, 0, 0.108814, 0, 0.0363762, 0, 0.0217084, 0, 0.0217784, 0, 0.0379644, 0, 0.0342662, 0, 0.0311304, 0, 0.0401337, 0, 0.0419444, 0, 0.0527232, 0, 0.0427343, 0, 0.0539223, 0, 0.0679392, 0.0791275, 0.0887525, 0.0778766, 0.0845525, 0.0826954, 0.065352, 0.0929196, 0.097882, 0.0824936, 0.0832312, 0.109067, 0.170409, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. };
//double MaxCluCoeff_26l[52] = { 0, 0.0877179, 0, 0.108814, 0, 0.0363762, 0, 0.0217084, 0, 0.0217784, 0, 0.0379644, 0, 0.0342662, 0, 0.0311304, 0, 0.0401337, 0, 0.0419444, 0, 0.0527232, 0, 0.0427343, 0, 0.0539223, 0, 0.0679392, 0.0791275, 0.0887525, 0.0778766, 0.0845525, 0.0826954, 0.065352, 0.0929196, 0.097882, 0.0824936, 0.0832312, 0.109067, 0.170409, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. };
//coefficient from C3d pions
double C3dCluCoeff_26l[52] = { 0, 0.0496421, 0, 0.0406136, 0, 0.0841348, 0, 0.0510543, 0, 0.0733869, 0, 0.0307511, 0, 0.0641894, 0, 0.0855152, 0, 0.0330387, 0, 0.0623249, 0, 0.0673885, 0, 0.067226, 0, 0.0921128, 0, 0.0835333, 0.103395, 0.0981839, 0.104409, 0.103367, 0.113767, 0.0369492, 0.0989041, 0.127434, 0.143822, 0.0649015, 0.0555431, 0.18853, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. };
//coefficient from photons - most energetic 2Dcluster/layer - 1.6<|eta|<2.8 - Seed Threshold = 10 mips, TrgHit Threshold = 5 mips
double MaxCluCoeffGamma_26l[52] = { 0, 0.0738292, 0, 0.0317296, 0, 0.0205507, 0, 0.0251915, 0, 0.0201699, 0, 0.0211367, 0, 0.0235583, 0, 0.0240063, 0, 0.0254955, 0, 0.0327115, 0, 0.0368678, 0, 0.040335, 0, 0.0354273, 0, 0.0667381, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. , 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
double MaxCluCoeffGamma_26l[52] = { 0., -0.0219692, 0., 0.060302, 0., 0.0294461, 0., 0.0273296, 0., 0.0235929, 0., 0.0276696, 0., 0.0271923, 0., 0.0291319, 0., 0.0331845, 0., 0.0338663, 0., 0.0441056, 0., 0.0134572, 0., 0.00226337, 0., 0.080997, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. , 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
for(int i=0; i<40; i++){
// %%%% For pions
// E += dep[i] * MaxCluCoeff_26l[i];
// E += dep[i] * C3dCluCoeff_26l[i];
// %%%% For photons
E += dep[i] * MaxCluCoeffGamma_26l[i];
//std::cout << " dep[i]*a[i] = " << dep[i] <<" * " << MaxCluCoeff_26l[i] << " = " <<E << std::endl;
std::cout << " dep[i]*a[i] = " << dep[i] <<" * " << C3dCluCoeff_26l[i] << " = " <<E << std::endl;
}
return E;
}
src/trgAnalysis.cpp
View file @ e91c0f5f
......@@ -42,6 +42,7 @@
#include "AnalysisOccupancy.hh"
#include "CoefficientDetTrg.hh"
#include "CoefficientDetClu.hh"
#include "CoefficientDetClu3D.hh"
#include "ReadEventTruthFile.hh"
#include "ReadEventTruthPu.hh"
......@@ -121,6 +122,7 @@ public:
std::cout << " 1024 - GammaShowerInfo" << std::endl;
std::cout << " 2048 - CoefficientDetTrg" << std::endl;
std::cout << " 4096 - CoefficientDetClu" << std::endl;
std::cout << " 8192 - CoefficientDetClu3D" << std::endl;
std::cout << "-b <evtBegin> First event to analyse, default " << evtBegin << std::endl;
std::cout << "-d <nDebug> Controls debug printout level, default " << nDebug << std::endl;
......@@ -354,6 +356,7 @@ int main(int argc, char *const *argv) {
if((nAnalysis&(1<< 10))!=0) vAnalysis.push_back(new GammaShowerInfo(sRoot,sOut));
if((nAnalysis&(1<< 11))!=0) vAnalysis.push_back(new CoefficientDetTrg(sRoot,sOut));
if((nAnalysis&(1<< 12))!=0) vAnalysis.push_back(new CoefficientDetClu(sRoot,sOut));
if((nAnalysis&(1<< 13))!=0) vAnalysis.push_back(new CoefficientDetClu3D(sRoot,sOut));
for(unsigned a(0);a<vAnalysis.size();a++) {
......
Markdown is supported
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