Commit fd5c2663 authored by dauncey's avatar dauncey

Trying to get update

parents afedd59b 42dce392
......@@ -105,8 +105,9 @@ public:
return s;
}
virtual void write(unsigned evt){;}
//virtual void write(unsigned evt){;}
//
protected:
const double fTwoPi;
......
#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
This diff is collapsed.
#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());