Skip to content

Commit

Permalink
Update the macro codes for IsoTrack and Muon analysis
Browse files Browse the repository at this point in the history
  • Loading branch information
Sunanda committed May 21, 2018
1 parent 7ba18fe commit dfe8038
Show file tree
Hide file tree
Showing 4 changed files with 495 additions and 197 deletions.
168 changes: 125 additions & 43 deletions Calibration/HcalCalibAlgos/macros/AnalyzeLepTree.C
View file
Expand Up @@ -13,14 +13,15 @@
// 2 : 0 ignores the information on track momentum
// 1 separate plots for certain momentum range
// (the range depends on ieta)
// 3 : 0 separate plot for each depth
// 3-4: 0 separate plot for each depth
// 1 sums up all depths
// 4-5: 0 no check on iphi
// 2 collapse depth
// 5-6: 0 no check on iphi
// 1 separate plot for each iphi
// 2 separate plot for each RBX
// 3 exclude the RBX specified by bits 6-10
// 6-10: RBX # to be excluded (maximum 5 bits needed for RBX)
// 11: 0 varying ranges of p depending on |ieta|
// 7-11: RBX # to be excluded (maximum 5 bits needed for RBX)
// 12: 0 varying ranges of p depending on |ieta|
// 1 constant p ranges
// modeLHC (integer) specifies the detector condition
// 0 Run1 detector (till 2016)
Expand Down Expand Up @@ -95,12 +96,13 @@ private:
TCanvas* plotHisto(TH1D* hist);
void plot2DHisto(std::map<unsigned int,TH2D*> hists,
std::vector<TCanvas*>& cvs, bool save);
int getCollapsedDepth(int eta, int phi, int depth);
int getRBX(int eta);
int getPBin(int eta);
int getVxBin();
int getDepthBin(int depth);
int getPhiBin(int eta);
int nDepthBins(int eta, int phi);
int nDepthBins(int eta, int phi, int modeLHC);
int nPhiBins(int eta);
int nPBins(int eta);
int nVxBins();
Expand Down Expand Up @@ -139,7 +141,7 @@ private:

static const int etamax_=26, npbin_=9, nvbin_=6;
static const bool debug_=false;
int mode_, modeLHC_, exRBX_, kphi_;
int mode_, modeLHC_, exRBX_, kphi_, kdepth_;
std::vector<int> npvbin_, iprange_;
std::vector<double> prange_[5];
std::map<unsigned int, TH1D*> h_p_, h_nv_;
Expand Down Expand Up @@ -195,8 +197,9 @@ void AnalyzeLepTree::Init(TTree *tree) {
// Init() will be called many times when running on PROOF
// (once per file to be processed).

exRBX_ = (mode_/64)%32;
kphi_ = (mode_/16)%4;
exRBX_ = (mode_/128)%32;
kphi_ = (mode_/32)%4;
kdepth_ = (mode_/8)%4;
if ((mode_%2) == 0) std::cout << "Produce set of histograms of energy, "
<< " active length, active length corrected "
<< "energy for 3 types" << std::endl;
Expand All @@ -210,9 +213,10 @@ void AnalyzeLepTree::Init(TTree *tree) {
<< "momentum" << std::endl;
else std::cout << "Separate plots for certain momentum "
<< "range (the range depends on ieta)\n";
if (((mode_/8)%2) == 0) std::cout << "Generate separate plot for each depth"
if (kdepth_ == 0) std::cout << "Generate separate plot for each depth"
<< std::endl;
else std::cout << "Sums up all depths for plots\n";
else if (kdepth_ == 1) std::cout << "Sums up all depths for plots\n";
else std::cout << "Collapse depths to Run 1 scenario\n";
if (kphi_ == 0) std::cout << "Make no check on iphi" << std::endl;
else if (kphi_ == 1) std::cout << "Make separate plot for each iphi\n";
else if (kphi_ == 2) std::cout << "Make separate plot for each RBX\n";
Expand Down Expand Up @@ -301,11 +305,13 @@ void AnalyzeLepTree::Loop(Long64_t nmax, bool debug) {
Long64_t nentries = fChain->GetEntriesFast();
std::cout << "Number of entries: " << nentries << ":" << nmax << std::endl;
if (nmax > 0 && nmax < nentries) nentries = nmax;
const double ethr = 0.00001; // Threshold of 10 keV

Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries; jentry++) {
Long64_t ientry = LoadTree(jentry);
if (debug) std::cout << "Entry " << jentry << ":" << ientry << std::endl;
if ((jentry%1000000 == 0) || debug)
std::cout << "Entry " << jentry << ":" << ientry << std::endl;
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
// if (Cut(ientry) < 0) continue;
Expand All @@ -330,15 +336,15 @@ void AnalyzeLepTree::Loop(Long64_t nmax, bool debug) {
if (it5 != h_nv2_.end()) (it5->second)->Fill(t_nvtx);
} else {
if (phi > 0 && pbin >= 0 && vbin >= 0) {
if ((mode_/8)%2 == 0) {
if (kdepth_ == 0) {
for (unsigned int k=0; k<t_depth->size(); ++k) {
int depth = (*t_depth)[k];
unsigned int id = packID(zside,eta,phi,depth+1,vbin,pbin);
double ene = (*t_ene)[k];
double enec = (*t_enec)[k];
double charge = (*t_charge)[k];
double actl = (*t_actln)[k];
if (ene > 0 && actl > 0 && charge > 0) {
if (ene > ethr && actl > 0 && charge > 0) {
std::map<unsigned int,TH1D*>::iterator it1 = h_Energy_.find(id);
if (it1 != h_Energy_.end()) (it1->second)->Fill(ene);
std::map<unsigned int,TH1D*>::iterator it2 = h_Ecorr_.find(id);
Expand Down Expand Up @@ -367,7 +373,7 @@ void AnalyzeLepTree::Loop(Long64_t nmax, bool debug) {
}
}
}
} else {
} else if (kdepth_ == 1) {
double ene(0), enec(0), actl(0), charge(0);
unsigned int id = packID(zside,eta,phi,1,vbin,pbin);
for (unsigned int k=0; k<t_depth->size(); ++k) {
Expand All @@ -378,18 +384,52 @@ void AnalyzeLepTree::Loop(Long64_t nmax, bool debug) {
actl += (*t_actln)[k];
}
}
std::map<unsigned int,TH1D*>::iterator it1 = h_Energy_.find(id);
if (it1 != h_Energy_.end()) (it1->second)->Fill(ene);
std::map<unsigned int,TH1D*>::iterator it2 = h_Ecorr_.find(id);
if (it2 != h_Ecorr_.end()) (it2->second)->Fill(ene/actl);
std::map<unsigned int,TH1D*>::iterator it3 = h_EnergyC_.find(id);
if (it3 != h_EnergyC_.end()) (it3->second)->Fill(enec);
std::map<unsigned int,TH1D*>::iterator it4 = h_EcorrC_.find(id);
if (it4 != h_EcorrC_.end()) (it4->second)->Fill(enec/actl);
std::map<unsigned int,TH1D*>::iterator it5 = h_Charge_.find(id);
if (it5 != h_Charge_.end()) (it5->second)->Fill(charge);
std::map<unsigned int,TH1D*>::iterator it6 = h_Chcorr_.find(id);
if (it6 != h_Chcorr_.end()) (it6->second)->Fill(charge/actl);
if (ene > ethr && actl > 0 && charge > 0) {
std::map<unsigned int,TH1D*>::iterator it1 = h_Energy_.find(id);
if (it1 != h_Energy_.end()) (it1->second)->Fill(ene);
std::map<unsigned int,TH1D*>::iterator it2 = h_Ecorr_.find(id);
if (it2 != h_Ecorr_.end()) (it2->second)->Fill(ene/actl);
std::map<unsigned int,TH1D*>::iterator it3 = h_EnergyC_.find(id);
if (it3 != h_EnergyC_.end()) (it3->second)->Fill(enec);
std::map<unsigned int,TH1D*>::iterator it4 = h_EcorrC_.find(id);
if (it4 != h_EcorrC_.end()) (it4->second)->Fill(enec/actl);
std::map<unsigned int,TH1D*>::iterator it5 = h_Charge_.find(id);
if (it5 != h_Charge_.end()) (it5->second)->Fill(charge);
std::map<unsigned int,TH1D*>::iterator it6 = h_Chcorr_.find(id);
if (it6 != h_Chcorr_.end()) (it6->second)->Fill(charge/actl);
}
} else {
double ene[3], enec[3], actl[3], charge[3];
for (unsigned int k=0; k<3; ++k) {
ene[k] = enec[k] = actl[k] = charge[k] = 0;
}
for (unsigned int k=0; k<t_depth->size(); ++k) {
if ((*t_ene)[k] > 0 && (*t_actln)[k] > 0) {
int dep = (*t_depth)[k];
int depth = getCollapsedDepth(zside*eta, phi, dep) - 1;
ene[depth] += (*t_ene)[k];
enec[depth] += (*t_enec)[k];
charge[depth] += (*t_charge)[k];
actl[depth] += (*t_actln)[k];
}
}
for (int k=0; k<nDepthBins(eta,phi,0); ++k) {
if (ene[k] > ethr && actl[k] > 0 && charge[k] > 0) {
unsigned int id = packID(zside,eta,phi,k+1,vbin,pbin);
std::map<unsigned int,TH1D*>::iterator it1 = h_Energy_.find(id);
if (it1 != h_Energy_.end()) (it1->second)->Fill(ene[k]);
std::map<unsigned int,TH1D*>::iterator it2 = h_Ecorr_.find(id);
if (it2 != h_Ecorr_.end()) (it2->second)->Fill(ene[k]/actl[k]);
std::map<unsigned int,TH1D*>::iterator it3 = h_EnergyC_.find(id);
if (it3 != h_EnergyC_.end()) (it3->second)->Fill(enec[k]);
std::map<unsigned int,TH1D*>::iterator it4 = h_EcorrC_.find(id);
if (it4 != h_EcorrC_.end()) (it4->second)->Fill(enec[k]/actl[k]);
std::map<unsigned int,TH1D*>::iterator it5 = h_Charge_.find(id);
if (it5 != h_Charge_.end()) (it5->second)->Fill(charge[k]);
std::map<unsigned int,TH1D*>::iterator it6 = h_Chcorr_.find(id);
if (it6 != h_Chcorr_.end()) (it6->second)->Fill(charge[k]/actl[k]);
}
}
}
}
}
Expand All @@ -411,7 +451,7 @@ void AnalyzeLepTree::bookHisto() {
double prange4[npbin_] = {0,125,150,175,200,250,300,400,500};
double prangeX[npbin_] = {125,150,200,250,300,350,400,500};
for (int i=0; i<npbin_; ++i) {
if ((mode_/2048)%2 == 0) {
if ((mode_/4096)%2 == 0) {
prange_[0].push_back(prange0[i]);
prange_[1].push_back(prange1[i]);
prange_[2].push_back(prange2[i]);
Expand All @@ -432,7 +472,8 @@ void AnalyzeLepTree::bookHisto() {
for (int ieta=-etamax_; ieta<=etamax_; ++ieta) {
int eta = (ieta>0) ? ieta : -ieta;
if (eta != 0) {
int ndepth = ((mode_/8)%2 == 0) ? nDepthBins(eta, 63) : 1;
int ndepth = ((kdepth_ == 0) ? nDepthBins(eta, 63, modeLHC_) :
(kdepth_ == 1) ? 1 : nDepthBins(eta, 63, 0));
std::cout << "Eta " << ieta << " with " << nPhiBins(eta)
<< " phi bins " << ndepth << " maximum depths and "
<< nPBins(eta) << " p bins" << std::endl;
Expand Down Expand Up @@ -520,10 +561,12 @@ void AnalyzeLepTree::bookHisto() {
} else {
sprintf (phis, "All i#phi");
}
int ndepth = ((mode_/8)%2 == 0) ? nDepthBins(eta, phi0) : 1;
int ndepth = ((kdepth_ == 0) ? nDepthBins(eta, phi0, modeLHC_) :
(kdepth_ == 1) ? 1 : nDepthBins(eta, phi0, 0));
for (int depth=0; depth<ndepth; ++depth) {
char deps[20];
if ((mode_/8)%2 == 0) sprintf (deps, "Depth=%d", depth+1);
if (kdepth_ == 1) sprintf (deps, "all depths");
else sprintf (deps, "Depth=%d", depth+1);
for (int pbin=0; pbin<nPBins(eta); ++pbin) {
char ps[30];
if ((mode_/4)%2 == 1) {
Expand All @@ -544,32 +587,32 @@ void AnalyzeLepTree::bookHisto() {
}
unsigned int id = packID(zside,eta,phi,depth+1,vbin,pbin);
char name[100], title[200];
sprintf (name,"EdepE%dF%dD%dV%dP%d",ieta,phi,depth,vbin,pbin);
sprintf (name,"EdepE%dF%dD%dV%dP%d",ieta,phi,deps,vbin,pbin);
sprintf (title,"Deposited energy for %s %s %s %s %s (GeV)", etas, phis, deps, ps, vtx);
getBins(0, ieta, phi0, depth+1, nbin, xmax);
h_Energy_[id] = new TH1D(name,title,nbin,0.0,xmax);
++book1;
sprintf (name,"EcorE%dF%dD%dV%dP%d",ieta,phi,depth,vbin,pbin);
sprintf (name,"EcorE%dF%dD%dV%dP%d",ieta,phi,deps,vbin,pbin);
sprintf (title,"Active length corrected energy for %s %s %s %s %s (GeV/cm)", etas, phis, deps, ps, vtx);
getBins(1, ieta, phi0, depth+1, nbin, xmax);
h_Ecorr_[id] = new TH1D(name,title,nbin,0.0,xmax);
++book1;
sprintf (name,"EdepCE%dF%dD%dV%dP%d",ieta,phi,depth,vbin,pbin);
sprintf (name,"EdepCE%dF%dD%dV%dP%d",ieta,phi,deps,vbin,pbin);
sprintf (title,"Response Corrected deposited energy for %s %s %s %s %s (GeV)", etas, phis, deps, ps, vtx);
getBins(2, ieta, phi0, depth+1, nbin, xmax);
h_EnergyC_[id] = new TH1D(name,title,nbin,0.0,xmax);
++book1;
sprintf (name,"EcorCE%dF%dD%dV%dP%d",ieta,phi,depth,vbin,pbin);
sprintf (name,"EcorCE%dF%dD%dV%dP%d",ieta,phi,deps,vbin,pbin);
sprintf (title,"Response and active length corrected energy for %s %s %s %s %s (GeV/cm)", etas, phis, deps, ps, vtx);
getBins(3, ieta, phi0, depth+1, nbin, xmax);
h_EcorrC_[id] = new TH1D(name,title,nbin,0.0,xmax);
++book1;
sprintf (name,"ChrgE%dF%dD%dV%dP%d",ieta,phi,depth,vbin,pbin);
sprintf (name,"ChrgE%dF%dD%dV%dP%d",ieta,phi,deps,vbin,pbin);
sprintf (title,"Measured charge for %s %s %s %s %s (cm)", etas, phis, deps, ps, vtx);
getBins(4, ieta, phi0, depth+1, nbin, xmax);
h_Charge_[id] = new TH1D(name,title,nbin,0.0,xmax);
++book1;
sprintf (name,"ChcorE%dF%dD%dV%dP%d",ieta,phi,depth,vbin,pbin);
sprintf (name,"ChcorE%dF%dD%dV%dP%d",ieta,phi,deps,vbin,pbin);
sprintf (title,"Active length corrected charge for %s %s %s %s %s (cm)", etas, phis, deps, ps, vtx);
getBins(5, ieta, phi0, depth+1, nbin, xmax);
h_Chcorr_[id] = new TH1D(name,title,nbin,0.0,xmax);
Expand Down Expand Up @@ -616,7 +659,8 @@ void AnalyzeLepTree::writeHisto(std::string fname) {
phi0 = 4*iphi + 1;
phi = iphi + 1;
};
int ndepth = ((mode_/8)%2 == 0) ? nDepthBins(eta, phi0) : 1;
int ndepth = ((kdepth_ == 0) ? nDepthBins(eta, phi0, modeLHC_) :
(kdepth_ == 1) ? 1 : nDepthBins(eta, phi0, 0));
for (int depth=0; depth<ndepth; ++depth) {
for (int pbin=0; pbin<nPBins(eta); ++pbin) {
for (int vbin=0; vbin<nVxBins(); ++vbin) {
Expand Down Expand Up @@ -779,6 +823,44 @@ void AnalyzeLepTree::plot2DHisto(std::map<unsigned int,TH2D*> hists,
}
}

int AnalyzeLepTree::getCollapsedDepth(int eta, int phi, int dep) {
int depth = dep+1;
if (std::abs(eta) <= 14 || std::abs(eta) == 17) {
depth = 1;
} else if (std::abs(eta) == 15) {
if (modeLHC_ > 3) {
if (dep > 3) depth = 2;
else depth = 1;
}
} else if (std::abs(eta) == 16) {
if (modeLHC_ == 0 || (modeLHC_ == 3 && (phi < 63 || phi > 66 || eta < 0))) {
} else {
if (dep == 3) depth = 2;
else if (dep == 4) depth = 3;
}
} else if (std::abs(eta) < 26) {
if (modeLHC_ == 0 || (modeLHC_ == 3 && (phi < 63 || phi > 66 || eta < 0))) {
} else {
if (dep < 3) depth = 1;
else depth = 2;
}
} else if (std::abs(eta) == 26) {
if (modeLHC_ == 0 || (modeLHC_ == 3 && (phi < 63 || phi > 66 || eta < 0))) {
} else {
if (dep < 4) depth = 1;
else depth = 2;
}
} else {
if (modeLHC_ == 0 || (modeLHC_ == 3 && (phi < 63 || phi > 66 || eta < 0))) {
} else {
if (dep < 3) depth = 1;
else if (dep == 3) depth = 2;
else depth = 3;
}
}
return depth;
}

int AnalyzeLepTree::getRBX(int eta) {
int rbx(1);
int phi = (eta > 20) ? (2*t_iphi + 1) : (t_iphi + 1);
Expand Down Expand Up @@ -810,7 +892,7 @@ int AnalyzeLepTree::getVxBin() {
}

int AnalyzeLepTree::getDepthBin(int depth) {
int bin = ((mode_/8)%2 == 0) ? depth : 1;
int bin = (kdepth_ == 0) ? depth : 1;
return bin;
}

Expand All @@ -826,7 +908,7 @@ int AnalyzeLepTree::getPhiBin(int eta) {
return bin;
}

int AnalyzeLepTree::nDepthBins(int eta, int phi) {
int AnalyzeLepTree::nDepthBins(int eta, int phi, int modeLHC) {
// Run 1 scenario
int nDepthR1[29]={1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,3,1,2,2,2,2,2,2,2,2,2,3,3,2};
// Run 2 scenario from 2018
Expand All @@ -841,13 +923,13 @@ int AnalyzeLepTree::nDepthBins(int eta, int phi) {
// = 3 --> corresponds to 2017 (Plan 1)
// = 4 --> corresponds to Run4 (post LS3)
int nbin(0);
if (modeLHC_ == 0) {
if (modeLHC == 0) {
nbin = nDepthR1[eta-1];
} else if (modeLHC_ == 1) {
} else if (modeLHC == 1) {
nbin = nDepthR2[eta-1];
} else if (modeLHC_ == 2) {
} else if (modeLHC == 2) {
nbin = nDepthR3[eta-1];
} else if (modeLHC_ == 3) {
} else if (modeLHC == 3) {
if (phi > 0) {
if (eta >= 16 && phi >= 63 && phi <= 66) {
nbin = nDepthR2[eta-1];
Expand Down

0 comments on commit dfe8038

Please sign in to comment.