Skip to content

GitLab

Commit 3cb5dc8c authored by Luca Mastrolorenzo's avatar Luca Mastrolorenzo
Browse files
Options

1) Removed obsolete analysis (Pion and Tau are replaced by the... 

1) Removed obsolete analysis (Pion and Tau are replaced by the CoefficientDetClu3d in many aspects); 2) Updates on CoefficientDetClu3d and cluster calibration coefficient; 3) C2d Threhsold have been made input parameters in order to easily run on the batch system
parent 12ebe2b7
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 417 additions and 1594 deletions
interface/AnalysisPi.hh deleted 100644 → 0
View file @ 12ebe2b7
#ifndef AnalysisPi_HH
#define AnalysisPi_HH
#include <iostream>
#include <cassert>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
#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 AnalysisPi : public AnalysisBase {
public:
enum {
kCutBins=50
};
AnalysisPi(const std::string &sRoot="", const std::string &sOut="")
: AnalysisBase("AnalysisPi",sRoot,sOut) {
h_PiE = new TH1F("h_PiE","Pi Energy Distribution",100, 0,100);
h_PiPt = new TH1F("h_PiPt","Pi Pt Distribution",100, 0,100);
h_PiEta = new TH1F("h_PiEta","Pi #eta Distribution",100, -4, 4);
h_PiPhi = new TH1F("h_PiEPhi","Pi #phi Distribution",100, -4, 4);
h_DR = new TH1F("h_DR", "#DeltaR barycenter-hits", 200, 0, 50);
h_Rho = new TH1F("h_Rho", "#rho distance distribution",1000,0,100);
h_XY = new TH2F("h_XY", "x-y distance distribution",1000,-50,50,1000,-50,50);
h_SimHit_XY= new TH2F("h_SimHit_XY", "SimHit XY",3600,-180,180,3600,-180,180);
h_SimHit_XYscale = new TH2F("h_SimHit_XYscale", "SimHit XY scale",2000,-1,1,2000,-1,1);
h_dDaugSim_scale = new TH2F("h_dDaugSim_scale", "Projected footprint daughter-simHit",2000,-1,1,2000,-1,1);
//h_ClusterFootPrint = new TH2F("h_ClusterFootPrint", "Cluster footprint",2000,-20,20,2000,-20,20);
h_MaxSimHitE = new TH1F("h_MaxSimHitE", "masx sim hit energy in the cone", 1000, 0, 10);
h_zBarycenter = new TH1F("h_zBarycenter", "z position of the barycenter",200,300,500);
h_PiDecMode = new TH1F("h_PiDecMode", "Pi Hadronic decay mode", 6, 0, 6);
h_PiVisMass= new TH1F("h_PiVisMass", "Pi Hadronic visible mass", 100, 0, 2);
h_PiMass= new TH1F("h_PiMass", "Pi Lepton mass", 100, 0, 2);
h_PiInEndcap = new TH2F("h_PiInEndcap", "Number of taus in each endcap",3,0,3,10,0,10);
h_DiPiMass = new TH1F("h_DiPiMass", "di-#tau invariant mass distribution",120,0,120);
h_EtDensityLayer = new TH1F("h_EtDensityLayer", "Transverse energy density per layer",45,0,45);
h_Layer = new TH1F("h_Layer", "Involved Layer distribution",45,0,45);
h_HoE = new TH1F("h_HoE", "Hadronic/Electromagnetic energy distribution",200,0,1);
for(unsigned l(0);l<=Geometry::kNumberOfLayers;l++) {
hPiVsCone[l]=new TH2F((fName+sLayerLabel[l]+"DaugVsCone").c_str(),
(sLayerTitle[l]+";Cone deposit;Daughter deposit").c_str(),
100,0.0,500.0,100,0.0,500.0);
hDaugXY[l]=new TH2F((fName+sLayerLabel[l]+"DaugXY").c_str(),
(sLayerTitle[l]+";Daughter XY;Daughter XY").c_str(),
300,-150,150,300,-150,150);
hPiXY[l]=new TH2F((fName+sLayerLabel[l]+"PiXY").c_str(),
(sLayerTitle[l]+";Pi XY;Pi XY").c_str(),
300,-150,150,300,-150,150);
h_ClusterFootPrintEntries[l]=new TH2F((fName+sLayerLabel[l]+"FootPrintEntriesXY").c_str(),
(sLayerTitle[l]+";FootPrint Entries XY;FootPrint Entries XY").c_str(),
50,-20,20,50,-20,20);
h_ClusterFootPrint[l]=new TH2F((fName+sLayerLabel[l]+"FootPrintXY").c_str(),
(sLayerTitle[l]+";FootPrint XY;FootPrint XY").c_str(),
50,-20,20,50,-20,20);
h_ClusterFootPrintMIP[l]=new TH2F((fName+sLayerLabel[l]+"FootPrintMIPXY").c_str(),
(sLayerTitle[l]+";FootPrintMIP XY;FootPrintMIP XY").c_str(),
50,-20,20,50,-20,20);
h_ClusterFootPrintEtaPhiMIP[l]=new TH2F((fName+sLayerLabel[l]+"FootPrintEtaPhiMIP").c_str(),
(sLayerTitle[l]+";FootPrintEtaPhiMIP;FootPrintEtaPhiMIP").c_str(),
200,-1,1,200,-1,1);
h_ClusterFootPrintEtaPhiEntries[l]=new TH2F((fName+sLayerLabel[l]+"FootPrintEtaPhiEntries").c_str(),
(sLayerTitle[l]+";FootPrintEtaPhiEntries;FootPrintEtaPhiEntries").c_str(),
200,-1,1,200,-1,1);
h_ClusterFootPrintEtaPhiMIP_norm[l]=new TH2F((fName+sLayerLabel[l]+"FootPrintEtaPhi_norm").c_str(),
(sLayerTitle[l]+";FootPrintEtaPhi_norm;FootPrintEtaPhi_norm").c_str(),
200,-1,1,200,-1,1);
}
hPiXYtot=new TH2F("hPiXYtot", "Pi XY", 300,-150,150,300,-150,150);
hDaugXYtot=new TH2F("hDaugXYtot", "PiDaugters XY",300,-150,150,300,-150,150);
}
virtual ~AnalysisPi() {
}
virtual bool event(Event &event) {
std::cout << fName << " Event " << fEventNumber << std::endl;
// Loop over all cells
for(unsigned e(0);e<Geometry::kNumberOfEndcaps;e++) {
//for(unsigned e(0);e<1;e++) {
std::cout << "ENDCAP - " << e << std::endl;
std::vector<Particle> &vPi(event.pions(e));
std::cout << "e = " << e << " pion size " << vPi.size() <<std::endl;
double Et_inPiCone=0.;
for(unsigned i(0);i<vPi.size();i++) {
h_PiE->Fill(vPi[i].FourP().E());
h_PiPt->Fill(vPi[i].FourP().Pt());
h_PiEta->Fill(vPi[i].FourP().Eta());
h_PiPhi->Fill(vPi[i].FourP().Phi());
//std::cout << "pion " << i << " (E,pt,eta,phi) = " << vPi[i].FourP().E() << ", " << vPi[i].FourP().Pt() << ", " << vPi[i].FourP().Eta() << ", " << vPi[i].FourP().Phi() << std::endl;
}
h_PiInEndcap->Fill( e, vPi.size() );
double EmEt=0.;
double HadEt=0.;
for(unsigned l(0);l<Geometry::kNumberOfLayers;l++) {
// std::cout << "LAYER " << l << std::endl;
// Find projection of tau into layer
std::vector<Particle> vPiProj;
bool useful(false);
for(unsigned i(0);i<vPi.size();i++){
//std::cout << "vtauDecMode[i] = " << vPiDecMode[i] << "\n";
vPiProj.push_back(vPi[i].project(Geometry::layerZ(e,l)));
}
for(unsigned i(0);i<vPiProj.size();i++){
// std::cout << " loop proj = " << vPiProj[i].position().rho() << std::endl;
if( !(vPiProj.size()>0 && vPiProj[i].FourP().Pt()>0.0 && vPiProj[i].position().rho()<150.0 ) ) continue;
useful=true;
hPiXY[l]->Fill(vPiProj[i].position().x(), vPiProj[i].position().y());
hPiXYtot->Fill(vPiProj[i].position().x(), vPiProj[i].position().y());
}
if(useful){
double MIP_TOT_layer = 0.;
for(unsigned w(0);w<Geometry::numberOfWafers(l);w++) {
if(Geometry::validWafer(l,w)) {
for(unsigned c(0);c<Geometry::numberOfCells(l,w);c++) {
Point p(Geometry::waferCellPoint(e,l,w,c));
std::vector<SimHit> &vSimHit(event.simHits(e,l,w,c));
double sumCell(0.0);
for(unsigned i(0);i<vSimHit.size();i++){
sumCell+=vSimHit[i].mips();
}
MIP_TOT_layer+=sumCell;
}
}
}
// Energy sums
std::vector<double> sumCone(vPi.size(),0.0);
std::vector<double> sumDaug(vPi.size(),0.0);
for(unsigned w(0);w<Geometry::numberOfWafers(l);w++) {
if(Geometry::validWafer(l,w)) {
for(unsigned c(0);c<Geometry::numberOfCells(l,w);c++) {
// Position of sim hit
Point p(Geometry::waferCellPoint(e,l,w,c));
//std::cout << " coordinates of the point = "<< p.x()<< ", "<< p.y() << std::endl;
// Loop over simHits
std::vector<SimHit> &vSimHit(event.simHits(e,l,w,c));
double sumCell(0.0);
double sumCellEt(0.0);
//std::cout << " size GIUSTA " << vSimHit.size() << std::endl;
for(unsigned i(0);i<vSimHit.size();i++){
sumCell+=vSimHit[i].mips();
sumCellEt+=vSimHit[i].transverseEnergy();
//sumCellEt+=vSimHit[i].mips();
}
bool inDR = false;
if(sumCell>0){
for(unsigned ii(0);ii<vPiProj.size();ii++) {
TVector3 Barycenter(vPiProj[ii].position().x(), vPiProj[ii].position().y(), vPiProj[ii].position().z() );
Point pBarycenter(Barycenter.x(), Barycenter.y(), Barycenter.z() );
std::cout << "barycenter (x,y,z) " << pBarycenter.x() << ", " << pBarycenter.y() << ", " << pBarycenter.z() << std::endl;
h_zBarycenter->Fill(fabs(Barycenter.z()));
Point2D pBarycenterScale(pBarycenter.scale());
//double EmaxSimHit = 0.;
//int maxSimHitidx = 0;
std::cout << "simHit Size--->" << vSimHit.size() << std::endl;
if(vSimHit.size()>0){
double DR = TMath::Sqrt( ( pBarycenter.x()-p.x() )*( pBarycenter.x()-p.x() ) + ( pBarycenter.y()-p.y() ) * ( pBarycenter.y()-p.y() ) );
h_DR->Fill(DR);
std::cout << "barycenter eta = " << pBarycenter.eta() << " hit eta = " << p.eta() << std::endl;
std::cout << "barycenter phi = " << pBarycenter.phi() << " hit phi = " << p.phi() << std::endl;
if(DR<20){
inDR = true;
Point pCentredInMax(p-pBarycenter);
std::cout << "Max " << pBarycenter.x()<< " sim " << p.x() << " diff " << pCentredInMax.x() << std::endl;
std::cout << "endcap " << e << " pi " << ii << " layer = " << l << " simHit(x,y) position relative to the max = " << pCentredInMax.x() << ", " << pCentredInMax.y() << std::endl;
h_ClusterFootPrintEntries[Geometry::kNumberOfLayers]->Fill( pCentredInMax.x(), pCentredInMax.y() );
h_ClusterFootPrint[Geometry::kNumberOfLayers]->Fill( pCentredInMax.x(), pCentredInMax.y(), sumCellEt );
h_ClusterFootPrintEntries[l]->Fill( pCentredInMax.x(), pCentredInMax.y());
h_ClusterFootPrint[l]->Fill( pCentredInMax.x(), pCentredInMax.y(), sumCellEt );
h_ClusterFootPrintMIP[l]->Fill( pCentredInMax.x(), pCentredInMax.y(), sumCell );
h_ClusterFootPrintEtaPhiMIP[l]->Fill( p.eta()-pBarycenter.eta(), deltaPhi(p.phi(),pBarycenter.phi()), sumCell );
h_ClusterFootPrintEtaPhiEntries[l]->Fill( p.eta()-pBarycenter.eta(), deltaPhi(p.phi(),pBarycenter.phi()) );
h_ClusterFootPrintEtaPhiMIP_norm[l]->SetBinContent(
h_ClusterFootPrintEtaPhiMIP_norm[l]->GetXaxis()->FindBin( p.eta()-pBarycenter.eta()),
h_ClusterFootPrintEtaPhiMIP_norm[l]->GetYaxis()->FindBin( deltaPhi(p.phi(),pBarycenter.phi())),
sumCell/MIP_TOT_layer);
h_ClusterFootPrintMIP[Geometry::kNumberOfLayers]->Fill( pCentredInMax.x(), pCentredInMax.y(), sumCell );
h_ClusterFootPrintEtaPhiMIP[Geometry::kNumberOfLayers]->Fill( p.eta()-pBarycenter.eta(), deltaPhi(p.phi(),pBarycenter.phi()), sumCell );
h_ClusterFootPrintEtaPhiMIP_norm[Geometry::kNumberOfLayers]->SetBinContent(
h_ClusterFootPrintEtaPhiMIP_norm[Geometry::kNumberOfLayers]->GetXaxis()->FindBin( p.eta()-pBarycenter.eta()),
h_ClusterFootPrintEtaPhiMIP_norm[Geometry::kNumberOfLayers]->GetYaxis()->FindBin( deltaPhi(p.phi(),pBarycenter.phi())),
sumCell/MIP_TOT_layer);
h_ClusterFootPrintEtaPhiEntries[Geometry::kNumberOfLayers]->Fill( p.eta()-pBarycenter.eta(), deltaPhi(p.phi(),pBarycenter.phi()) );
//std::cout << "z center = " <<-1*TMath::Log( TMath::Tan(TMath::ATan( pCentredInMax.rho()/p.z() )/2) ) << std::endl;
Et_inPiCone+=sumCellEt;
}
}
}
h_SimHit_XY->Fill(p.x(), p.y());
Point2D ps(p.scale());
h_SimHit_XYscale->Fill(ps.x(),ps.y());
if(inDR==true){
h_EtDensityLayer->Fill(l, sumCellEt);
h_Layer->Fill(l);
if(l<=28)EmEt+=sumCellEt;
else if(l>28)HadEt+=sumCellEt;
std::cout << "layer, E, H " << l << ", "<< EmEt << ", " << HadEt << std::endl;
}
for(unsigned j(0);j<vPiProj.size();j++) {
bool nearCone(false);
//bool nearDaug(false);
// Delta-R to tau projection
double dR( vPiProj[j].position().deltaR(p) );
if(dR<0.2) nearCone=true;
if(nearCone) sumCone[j]+=sumCell;
}
}
}
}
}
for(unsigned i(0);i<sumCone.size();i++) {
std::cout << "Endcap, layer = " << e << ", " << l
<< ", pion " << i << " cone = " << sumCone[i] << std::endl;
hPiVsCone[l]->Fill(sumCone[i],sumDaug[i]);
hPiVsCone[Geometry::kNumberOfLayers]->Fill(sumCone[i],sumDaug[i]);
}
}
}
std::cout << "-----------------> Pi Et estimated "<< Et_inPiCone << std::endl;
if(EmEt>0){
double ratioHoE = HadEt/EmEt;
std::cout <<"Final E, H" << EmEt<< ", "<< HadEt << ", "<< ratioHoE << std::endl;
h_HoE->Fill(ratioHoE);
}
}
return true;
}
private:
TH2F *hPiVsCone[Geometry::kNumberOfLayers+1];
TH2F *hDaugXY[Geometry::kNumberOfLayers+1];
TH2F *hPiXY[Geometry::kNumberOfLayers+1];
TH2F *h_ClusterFootPrintEntries[Geometry::kNumberOfLayers+1];
TH2F *h_ClusterFootPrint[Geometry::kNumberOfLayers+1];
TH2F *h_ClusterFootPrintMIP[Geometry::kNumberOfLayers+1];
TH2F *h_ClusterFootPrintEtaPhiMIP[Geometry::kNumberOfLayers+1];
TH2F *h_ClusterFootPrintEtaPhiEntries[Geometry::kNumberOfLayers+1];
TH2F *h_ClusterFootPrintEtaPhiMIP_norm[Geometry::kNumberOfLayers+1];
TH1F *h_EtDensityLayer;
TH1F *h_Layer;
TH1F *h_HoE;
TH2F *hDaugXYtot;
TH2F *hPiXYtot;
TH1F *h_PiE;
TH1F *h_PiPt;
TH1F *h_PiEta;
TH1F *h_PiPhi;
TH1F *h_DR;
TH1F *h_Rho;
TH2F *h_XY;
TH2F *h_SimHit_XY;
TH2F *h_SimHit_XYscale;
TH2F *h_dDaugSim_scale;
TH1F *h_zBarycenter;
TH1F *h_MaxSimHitE;
TH1F *h_PiDecMode;
TH1F *h_PiVisMass;
TH1F *h_PiMass;
TH2F *h_PiInEndcap;
TH1F *h_DiPiMass;
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
interface/AnalysisTau.hh deleted 100644 → 0
View file @ 12ebe2b7
This diff is collapsed.
interface/AnalysisTauReso.hh deleted 100644 → 0
View file @ 12ebe2b7
#ifndef AnalysisTauReso_HH
#define AnalysisTauReso_HH
#include <iostream>
#include <cassert>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
#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 AnalysisTauReso : public AnalysisBase {
public:
enum {
kCutBins=50
};
AnalysisTauReso(const std::string &sRoot="", const std::string &sOut="")
: AnalysisBase("AnalysisTauReso",sRoot,sOut) {
h_TauE = new TH1F("h_TauE","Tau Energy Distribution",100, 0,100);
h_TauPt = new TH1F("h_TauPt","Tau Pt Distribution",100, 0,100);
h_TauEta = new TH1F("h_TauEta","Tau #eta Distribution",100, -4, 4);
h_TauPhi = new TH1F("h_TauEPhi","Tau #phi Distribution",100, -4, 4);
h_TauRes_1p = new TH1F("h_TauRes_1p","Tau Energy Response - 1p",100, -1.5, 1.5);
h_TauVisRes_1p = new TH1F("h_TauVisRes_1p","Visible Tau Energy Response - 1p",100, -1.5, 1.5);
h_TauVisResRel_1p = new TH1F("h_TauVisResRel_1p","Visible Tau Relative Energy Response - 1p",100, -1.5, 1.5);
h_TauRes_1ppz = new TH1F("h_TauRes_1ppz","Tau Energy Response - 1ppz",100, -1.5, 1.5);
h_TauVisRes_1ppz = new TH1F("h_TauVisRes_1ppz","Visible Tau Energy Response - 1ppz",100, -1.5, 1.5);
h_TauVisResRel_1ppz = new TH1F("h_TauVisResRel_1ppz","Visible Tau Relative Energy Response - 1ppz",100, -1.5, 1.5);
h_TauRes_3p = new TH1F("h_TauRes_3p","Tau Energy Response - 3p",100, -1.5, 1.5);
h_TauVisRes_3p = new TH1F("h_TauVisRes_3p","Visible Tau Energy Response - 3p",100, -1.5, 1.5);
h_TauVisResRel_3p = new TH1F("h_TauVisResRel_3p","Visible Tau Relative Energy Response - 3p",100, -1.5, 1.5);
h_DR = new TH1F("h_DR", "delta R tau-daugthers", 200, 0, 50);
h_Rho = new TH1F("h_Rho", "#rho distance distribution",1000,0,100);
h_XY = new TH2F("h_XY", "x-y distance distribution",1000,-50,50,1000,-50,50);
h_SimHit_XY= new TH2F("h_SimHit_XY", "SimHit XY",3600,-180,180,3600,-180,180);
h_SimHit_XYscale = new TH2F("h_SimHit_XYscale", "SimHit XY scale",2000,-1,1,2000,-1,1);
h_dDaugSim_scale = new TH2F("h_dDaugSim_scale", "Projected footprint daughter-simHit",2000,-1,1,2000,-1,1);
//h_ClusterFootPrint = new TH2F("h_ClusterFootPrint", "Cluster footprint",2000,-20,20,2000,-20,20);
h_MaxSimHitE = new TH1F("h_MaxSimHitE", "masx sim hit energy in the cone", 1000, 0, 10);
h_zBarycenter = new TH1F("h_zBarycenter", "z position of the barycenter",200,300,500);
h_TauDecMode = new TH1F("h_TauDecMode", "Tau Hadronic decay mode", 6, 0, 6);
h_TauVisMass= new TH1F("h_TauVisMass", "Tau Hadronic visible mass", 100, 0, 2);
h_TauMass= new TH1F("h_TauMass", "Tau Lepton mass", 100, 0, 2);
h_TauInEndcap = new TH2F("h_TauInEndcap", "Number of taus in each endcap",3,0,3,10,0,10);
h_DiTauMass = new TH1F("h_DiTauMass", "di-#tau invariant mass distribution",120,0,120);
h_EtDensityLayer = new TH1F("h_EtDensityLayer", "Transverse energy density per layer",45,0,45);
h_Layer = new TH1F("h_Layer", "Involved Layer distribution",45,0,45);
h_HoE = new TH1F("h_HoE", "Hadronic/Electromagnetic energy distribution",200,0,1);
for(unsigned l(0);l<=Geometry::kNumberOfLayers;l++) {
hDaugVsCone[l]=new TH2F((fName+sLayerLabel[l]+"DaugVsCone").c_str(),
(sLayerTitle[l]+";Cone deposit;Daughter deposit").c_str(),
100,0.0,500.0,100,0.0,500.0);
hDaugXY[l]=new TH2F((fName+sLayerLabel[l]+"DaugXY").c_str(),
(sLayerTitle[l]+";Daughter XY;Daughter XY").c_str(),
300,-150,150,300,-150,150);
hTauXY[l]=new TH2F((fName+sLayerLabel[l]+"TauXY").c_str(),
(sLayerTitle[l]+";Tau XY;Tau XY").c_str(),
300,-150,150,300,-150,150);
h_ClusterFootPrintEntries[l]=new TH2F((fName+sLayerLabel[l]+"FootPrintEntries XY").c_str(),
(sLayerTitle[l]+";FootPrint Entries XY;FootPrint Entries XY").c_str(),
50,-20,20,50,-20,20);
h_ClusterFootPrint[l]=new TH2F((fName+sLayerLabel[l]+"FootPrint XY").c_str(),
(sLayerTitle[l]+";FootPrint XY;FootPrint XY").c_str(),
50,-20,20,50,-20,20);
}
hTauXYtot=new TH2F("hTauXYtot", "Tau XY", 300,-150,150,300,-150,150);
hDaugXYtot=new TH2F("hDaugXYtot", "TauDaugters XY",300,-150,150,300,-150,150);
}
virtual ~AnalysisTauReso() {
}
virtual bool event(Event &event) {
std::cout << fName << " Event " << fEventNumber << std::endl;
// Loop over all cells
for(unsigned e(0);e<Geometry::kNumberOfEndcaps;e++) {
std::vector<Particle> &vTau(event.taus(e));
std::vector< std::vector<Particle> > &vTauDaughter(event.tauDaughters(e));
std::vector< std::vector<Particle> > &vTauGrandChildren(event.tauGrandChildren(e));
std::vector< int > vTauDecMode;
std::vector< int > vTauVisEt;
// double Et_inTauCone=0.;
for(unsigned i(0);i<vTau.size();i++) {
h_TauE->Fill(vTau[i].FourP().E());
h_TauPt->Fill(vTau[i].FourP().Pt());
h_TauEta->Fill(vTau[i].FourP().Eta());
h_TauPhi->Fill(vTau[i].FourP().Phi());
int nPi=0;
int nPiZero=0;
int nRho=0;
int nRhoZero=0;
int nA1=0;
int TauDecMode = 0;
TLorentzVector Daug4P;
TLorentzVector VisDaug4P;
for(unsigned j(0);j<vTauDaughter[i].size();j++) {
std::cout << " -- daughter " << j << " pdgId: " << vTauDaughter[i][j].pdgId();
std::cout << " (E,pt,eta,phi) = " << vTauDaughter[i][j].FourP().E() << ", " << vTauDaughter[i][j].FourP().Pt();
std::cout << ", " << vTauDaughter[i][j].FourP().Eta() << ", " << vTauDaughter[i][j].FourP().Phi();
std::cout << " (Px,Py,Pz)" << vTauDaughter[i][j].FourP().Px() << ", " << vTauDaughter[i][j].FourP().Py() << ", " << vTauDaughter[i][j].FourP().Pz() << std::endl;
TLorentzVector tmpDaug4P;
tmpDaug4P.SetPxPyPzE(vTauDaughter[i][j].FourP().Px(), vTauDaughter[i][j].FourP().Py(), vTauDaughter[i][j].FourP().Pz(), vTauDaughter[i][j].FourP().E());
Daug4P += tmpDaug4P;
if( abs( vTauDaughter[i][j].pdgId() ) == 211 ){nPi = nPi+1; }
if( !( abs( vTauDaughter[i][j].pdgId() ) == 12 || abs( vTauDaughter[i][j].pdgId() ) == 14 || abs( vTauDaughter[i][j].pdgId() ) == 16 ) )VisDaug4P += tmpDaug4P;
}
for(unsigned r(0);r<vTauGrandChildren[i].size();r++) {
std::cout << " == grandchildren " << r << " pdgId: " << vTauGrandChildren[i][r].pdgId() << " (E,pt,eta,phi) = " << vTauGrandChildren[i][r].FourP().E() << ", " << vTauGrandChildren[i][r].FourP().Pt() << ", " << vTauGrandChildren[i][r].FourP().Eta() << ", " << vTauGrandChildren[i][r].FourP().Phi();
std::cout << " (Px,Py,Pz)" << vTauGrandChildren[i][r].FourP().Px() << ", " << vTauGrandChildren[i][r].FourP().Py() << ", " << vTauGrandChildren[i][r].FourP().Pz() << std::endl;
if( abs( vTauGrandChildren[i][r].pdgId() ) == 111 ){nPiZero = nPiZero+1.;}
if( abs( vTauGrandChildren[i][r].pdgId() ) == 213 ){nRho = nRho+1.; }
if( abs( vTauGrandChildren[i][r].pdgId() ) == 113 ){nRhoZero = nRhoZero+1.; }
if( abs( vTauGrandChildren[i][r].pdgId() ) == 20113 ){nA1 = nA1+1.; }
}
h_TauVisMass->Fill(VisDaug4P.M());
h_TauMass->Fill(Daug4P.M());
if(nPi==1 && nPiZero==0 && VisDaug4P.M() < 0.2){
TauDecMode = 0;
}
//---------- 1prong+1piZero -------1-prong+2piZero -
else if( VisDaug4P.M() > 0.2 && ( (nPi==1 && nPiZero==1) || (nRho==1 && nPiZero==1) || (nRho==1 && nPiZero==0 && nPi==0 && nA1==0) ) ){
TauDecMode = 1;
}
else if( VisDaug4P.M() > 0.4 && ( nRhoZero==1 && nPi>=1) ){
TauDecMode = 4;
}
else{ TauDecMode = 5;}
h_TauDecMode->Fill( TauDecMode );
vTauDecMode.push_back( TauDecMode );
vTauVisEt.push_back(VisDaug4P.Et());
}
for(unsigned t(0);t<vTau.size();t++) {
double Et_tot = 0.;
if(!(fabs(vTau[t].FourP().Eta())>1.47 && fabs(vTau[t].FourP().Eta())<3.0 && vTau[t].FourP().Pt()>40 )) continue;
for(unsigned l(0);l<Geometry::kNumberOfLayers;l++) {
//std::cout << "LAYER " << l << std::endl;
// Find projection of tau into layer
Particle TauProj = vTau[t].project(Geometry::layerZ(e,l));
double Et_tot_tmp=0.;
for(unsigned w(0);w<Geometry::numberOfWafers(l);w++) {
if(Geometry::validWafer(l,w)) {
for(unsigned c(0);c<Geometry::numberOfCells(l,w);c++) {
// Position of sim hit
Point p(Geometry::waferCellPoint(e,l,w,c));
//std::cout << " coordinates of the point = "<< p.x()<< ", "<< p.y() << std::endl;
// Loop over simHits
std::vector<SimHit> &vSimHit(event.simHits(e,l,w,c));
double sumCell(0.0);
double sumCellEt(0.0);
//std::cout << " size GIUSTA " << vSimHit.size() << std::endl;
for(unsigned i(0);i<vSimHit.size();i++){
sumCell+=vSimHit[i].mips();
sumCellEt+=vSimHit[i].transverseEnergy();
//sumCellEt+=vSimHit[i].mips();
}
//double DR = TMath::Sqrt( deltaPhi(p.phi(),TauProj.position().phi()) * deltaPhi(p.phi(),TauProj.position().phi()) + (p.eta() - TauProj.position().eta())*(p.eta() - TauProj.position().eta() ) );
double DR = TMath::Sqrt( (p.x() - TauProj.position().x())*(p.x() - TauProj.position().x() ) + (p.y() - TauProj.position().y())*(p.y() - TauProj.position().y() ) );
if(sumCellEt>0 && DR<20 ){
Et_tot_tmp+=sumCellEt;
}
}
}
}
Et_tot+=Et_tot_tmp;
}//loop over tau projection x layer
std::cout << "Eta " << vTau[t].FourP().Eta() << " Tau Et = "<< vTau[t].FourP().Et()<< " tau estimated energy = " << Et_tot << std::endl;
std::cout << " Visible Tau Et = "<< vTauVisEt[t] << " tau estimated energy = " << Et_tot << std::endl;
if(vTauDecMode[t]==0){
h_TauRes_1p->Fill(( vTau[t].FourP().Et()-Et_tot) /vTau[t].FourP().Et());
h_TauVisRes_1p->Fill( ( vTauVisEt[t]-Et_tot ) / vTauVisEt[t] );
h_TauVisResRel_1p->Fill( Et_tot / vTauVisEt[t] );
}
if(vTauDecMode[t]==1){
h_TauRes_1ppz->Fill(( vTau[t].FourP().Et()-Et_tot) /vTau[t].FourP().Et());
h_TauVisRes_1ppz->Fill( ( vTauVisEt[t]-Et_tot ) / vTauVisEt[t] );
h_TauVisResRel_1ppz->Fill( Et_tot / vTauVisEt[t] );
}
if(vTauDecMode[t]==4){
h_TauRes_3p->Fill(( vTau[t].FourP().Et()-Et_tot) /vTau[t].FourP().Et());
h_TauVisRes_3p->Fill( ( vTauVisEt[t]-Et_tot ) / vTauVisEt[t] );
h_TauVisResRel_3p->Fill( Et_tot / vTauVisEt[t] );
}
}
h_TauInEndcap->Fill( e, vTau.size() );
if(vTau.size()>1){
std::cout << "E1 = " << vTau[0].FourP().E() << " E2 = " << vTau[1].FourP().E() << "\n";
std::cout << "px1 = " << vTau[0].FourP().Px() << " px2 = " << vTau[1].FourP().Px() << "\n";
std::cout << "py1 = " << vTau[0].FourP().Py() << " py2 = " << vTau[1].FourP().Py() << "\n";
std::cout << "pz1 = " << vTau[0].FourP().Pz() << " pz2 = " << vTau[1].FourP().Pz() << "