//#define AnalysisBase_cxx
#include "AnalysisBase.h"
AnalysisBase::AnalysisBase(TTree *tree) : fChain(0)
{
// Notes
// For Rz alignment, one should SUBTRACT the slope found from the deltaX vs Y plot.
// For Z alignment, one should ADD the slope found from the deltaX vs Y plot.
// For Ry: If slope of DX vs X is negative (positive), reduce (increase) Ry
//
//these have to be able to change so we have to make sure there is only 1 copy of them, so global scope would be bad news
minDistFromHole = m_minDistFromHole;
holeSector = m_board.Contains("D") && (m_sector=="1" || m_sector=="2" || m_sector == "3");
if(!holeSector) minDistFromHole = -9999;
stripPitch = 0.190;
z_DUT = 370; //This needs to be set for different boards
Rz = -0.0218;
Ry = 0.00;
dxWin = 0.25;
xGloOff = -9.2;
yGloOff = -7.5;
xOff = -54.78;
correctForZRotation = true;
channelOffset = 0.0;
xLeftHole = 999;
xRightHole = -999;
stripGap[0] = 0;
stripGap[1] = 0;
stripGap[2] = 0;
stripGap[3] = 0;
nDeadRegion = 0;
if( m_board.Contains("D")) {
stripPitch = 0.095;
channelOffset = 512.0;
//stripGap[0]=19;
//stripGap[1]=39;
//stripGap[2]=58;
//stripGap[3]=0;
stripGap[0]=58;
stripGap[1]=39;
stripGap[2]=19;
stripGap[3]=0;
if(m_board.Contains("D5")) {
z_DUT = 1100.0;
//Rz = -0.0230;
}else if(m_board.Contains("D7")) {
z_DUT = 1123.0;
if(m_sector=="5") z_DUT = 1087.0;
if(m_sector=="1") z_DUT = 1040.0;
//Ry = -0.008;
}
}
yInt1[0] = 2.8; yInt1[1] = 3.4;
yInt2[0] = 3.4; yInt2[1] = 4.1;
yInt3[0] = 1.0; yInt3[1] = 2.3;
if(m_board.Contains("A8")){
z_DUT = 1100.0;
Rz = -0.0286;
yInt1[0] = 1.8; yInt1[1] = 2.0;
yInt2[0] = 2.0; yInt2[1] = 2.4;
yInt3[0] = 1.0; yInt3[1] = 2.3;
if(m_sector=="3" || m_sector=="6") {
channelOffset = -128.0*2+22;
}else{
channelOffset = -128.0;
}
if(m_scanType.Contains("Angle")){
if(m_angle == "10") {
Ry = 9.3*(TMath::Pi()/180.0);
z_DUT = z_DUT - 18.0;
} else if(m_angle == "20") {
Ry = 19*(TMath::Pi()/180.0);
z_DUT = z_DUT - 18.0;
}
}
}else if( m_board.Contains("A4")){
z_DUT = 381.0;
Rz = -0.0238;
if(m_sector=="3" || m_sector=="6") {
channelOffset = -128.0*2.0+19;
}else{
channelOffset = -128.0;
}
}else if( m_board.Contains("A6")){
yInt1[0] = 2.5; yInt1[1] = 3.1;
yInt2[0] = 3.1; yInt2[1] = 3.7;
yInt3[0] = 1.0; yInt3[1] = 2.3;
if(m_sector=="3" || m_sector=="6") {
channelOffset = -128.0*2+22;
}else{
channelOffset = -128.0;
}
}else if( m_board.Contains("A2") || m_board.Contains("A1")){
z_DUT = 310;
yInt1[0] = -5.0; yInt1[1] = -2.0;
yInt2[0] = -2.0; yInt2[1] = 1.0;
yInt3[0] = 1.0; yInt3[1] = 5.0;
}
// Board A1 alignment is really messed up for a few runs! ADHOC here!
if(m_board.Contains("A1")){// && (m_bias=="300" || m_bias=="350" || m_bias=="150" || m_bias==")){
dxWin = 1.0;
correctForZRotation = false;
}
holeQuadPar[0] = 0;
holeQuadPar[1] = 0;
holeQuadPar[2] = 0;
removeTracksInHole = removeTracksInHoleDef;
// Only look to remove tacks from hole area if D-type and in sectors 1, 2, or 3.
if(!m_board.Contains("D") || !(m_sector=="1" || m_sector=="2" || m_sector=="3") ) removeTracksInHole = false;
polarity = 1.0;
if(isPType) polarity = -1;
}
AnalysisBase::~AnalysisBase()
{
if (!fChain) return;
delete fChain->GetCurrentFile();
}
Int_t AnalysisBase::GetEntry(Long64_t entry)
{
// Read contents of entry.
if (!fChain) return 0;
return fChain->GetEntry(entry);
}
Long64_t AnalysisBase::LoadTree(Long64_t entry)
{
// Set the environment to read one entry
if (!fChain) return -5;
Long64_t centry = fChain->LoadTree(entry);
if (centry < 0) return centry;
if (fChain->GetTreeNumber() != fCurrent) {
fCurrent = fChain->GetTreeNumber();
Notify();
}
return centry;
}
void AnalysisBase::Init(TTree *tree)
{
// The Init() function is called when the selector needs to initialize
// a new tree or chain. Typically here the branch addresses and branch
// pointers of the tree will be set.
// It is normally not necessary to make changes to the generated
// code, but the routine can be extended by the user if needed.
// Init() will be called many times when running on PROOF
// (once per file to be processed).
// Set object pointer
vec_trk_x = 0;
vec_trk_y = 0;
vec_trk_tx = 0;
vec_trk_ty = 0;
vec_trk_chi2ndf = 0;
// Set branch addresses and branch pointers
if (!tree) return;
fChain = tree;
fCurrent = -1;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("clusterNumberPerEvent", &clusterNumberPerEvent, &b_clusterNumberPerEvent);
fChain->SetBranchAddress("clustersTDC", &clustersTDC, &b_clustersTDC);
fChain->SetBranchAddress("timestamps", ×tamps, &b_timestamps);
fChain->SetBranchAddress("clustersPosition", clustersPosition, &b_clustersPosition);
fChain->SetBranchAddress("clustersSeedPosition", clustersSeedPosition, &b_clustersSeedPosition);
fChain->SetBranchAddress("clustersCharge", clustersCharge, &b_clustersCharge);
fChain->SetBranchAddress("clustersSize", clustersSize, &b_clustersSize);
fChain->SetBranchAddress("clustersSeedCharge", clustersSeedCharge, &b_clustersSeedCharge);
fChain->SetBranchAddress("clustersCharge2StripLeft", clustersCharge2StripLeft, &b_clustersCharge2StripLeft);
fChain->SetBranchAddress("clustersCharge1StripLeft", clustersCharge1StripLeft, &b_clustersCharge1StripLeft);
fChain->SetBranchAddress("clustersCharge1StripRight", clustersCharge1StripRight, &b_clustersCharge1StripRight);
fChain->SetBranchAddress("clustersCharge2StripRight", clustersCharge2StripRight, &b_clustersCharge2StripRight);
fChain->SetBranchAddress("n_tp3_tracks", &n_tp3_tracks, &b_n_tp3_tracks);
fChain->SetBranchAddress("vec_trk_x", &vec_trk_x, &b_vec_trk_x);
fChain->SetBranchAddress("vec_trk_y", &vec_trk_y, &b_vec_trk_y);
fChain->SetBranchAddress("vec_trk_tx", &vec_trk_tx, &b_vec_trk_tx);
fChain->SetBranchAddress("vec_trk_ty", &vec_trk_ty, &b_vec_trk_ty);
fChain->SetBranchAddress("vec_trk_chi2ndf", &vec_trk_chi2ndf, &b_vec_trk_chi2ndf);
fChain->SetBranchAddress("dtime", &dtime, &b_dtime);
Notify();
}
Bool_t AnalysisBase::Notify()
{
// The Notify() function is called when a new file is opened. This
// can be either for a new TTree in a TChain or when when a new TTree
// is started when using PROOF. It is normally not necessary to make changes
// to the generated code, but the routine can be extended by the
// user if needed. The return value is currently not used.
return kTRUE;
}
void AnalysisBase::Show(Long64_t entry)
{
// Print contents of entry.
// If entry is not specified, print current entry
if (!fChain) return;
fChain->Show(entry);
}
/*
Int_t AnalysisBase::Cut(Long64_t entry)
{
// This function may be called from Loop.
// returns 1 if entry is accepted.
// returns -1 otherwise.
return 1;
}
*/
void AnalysisBase::getRange(TH1 *h, float &lo, float& hi, float thresh, int nSkipMax){
int nfail = 0;
int iMaxBin = h->GetMaximumBin();
lo = iMaxBin;
hi = iMaxBin;
double maxValue = h->GetBinContent(iMaxBin);
// Only for debugging....
//h->Draw();
//std::cout << "maxValue = " << iMaxBin << " " << maxValue << std::endl;
for(int j=iMaxBin;j>=0;j--){
double v = h->GetBinContent(j);
double r = v / maxValue;
//std::cout << "r = " << r << std::endl; // Only for debugging
if( r > thresh) {
lo = j ;
nfail = 0;
}else if( r <= thresh) {
nfail++;
}
if(nfail >= nSkipMax) break;
}
nfail = 0;
for(int j=iMaxBin;j<=h->GetNbinsX();j++){
double v = h->GetBinContent(j);
double r = v / maxValue;
//std::cout << "r = " << r << std::endl; // Only for debugging
if( r > thresh) {
hi = j;
nfail = 0;
}else if( r <= thresh) {
nfail++;
}
if(nfail >= nSkipMax) break;
}
lo = h->GetBinCenter(lo);
hi = h->GetBinCenter(hi);
//std::cout << "lo, hi = " << lo << " " << hi << std::endl; // only for debugging
return;
}
void AnalysisBase::getTDCBins(TProfile* h, float& lo, float& hi){
getRange(h,lo,hi,0.97,1);
lo = lo - 1.0 + 0.5;
hi = hi + 0.5;
return;
}
void AnalysisBase::getBeamLocation(TH1F *h, float &lo, float& hi){
//double dif = 10;
for(int j = 0; j<10; j++){
std::cout << "====> Finding Beam, iteration " << j+1 << std::endl;
getRange(h,lo,hi);
int dif = hi - lo + 1;
if(dif < 5){
int i1 = h->FindBin(lo);
int i2 = h->FindBin(hi);
for(int k=i1; k<=i2; k++){
h->SetBinContent(k,0);
}
}else{
break;
}
}
//h->Draw();
return;
}
double AnalysisBase::getXOffset(TH1F *h1w){
float xlo = 0, xhi=0;
getRange(h1w,xlo,xhi,0.1,1);
return (xlo+xhi)/2.0;
}
void AnalysisBase::getBeamLoc(){
TH1F* ha = new TH1F("ha","Strip # of cluster with track",512,0.0,512);
//TH1F* hb = new TH1F("hb","Strip # of cluster with track",512,0.0,512);
Long64_t nentries = fChain->GetEntriesFast();
float lo = 0, hi = 0;
Long64_t nbytes = 0, nb = 0;
std::cout << "======================================= " << std::endl;
std::cout << "getBeamLoc(): Determining Beam Position " << nentries << std::endl;
std::cout << "======================================= " << std::endl;
for (Long64_t jentry=0; jentry<max(50000,(int)nentries);jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
//cout << clusterNumberPerEvent << endl;
for(int j=0; j<min((int)clusterNumberPerEvent,10); j++){
//std::cout << clustersPosition[j] << std::endl;
//cout << polarity << " " << polarity*clustersCharge[j] << endl;
if(polarity*clustersCharge[j] < kClusterChargeMin) continue;
int iChan = clustersSeedPosition[j];
//if(j<500) cout << "iChan = " << iChan << clustersCharge[j] << " " << 5*noise[iChan] << endl;
if(polarity*clustersCharge[j]<4*noise[iChan]) continue;
if(clustersPosition[j]>0.1&&clustersSize[j]==1) ha->Fill(clustersPosition[j]);
if(clustersPosition[j]>0.1&&clustersSize[j]==2) ha->Fill(clustersPosition[j]);
}
}
int num = ha->GetEntries();
for(int i=0; i<1000; i++){
if(num < 1000) continue;
break;
}
if(num < 1000){
std::cout << "ERROR: Something wrong here, insufficient entries in GetBeamLoc(), nEntries = " << num << std::endl;
//exit(1);
}
//ha->Draw();
getBeamLocation(ha,lo,hi);
iLo = lo;
iHi = hi;
std::cout << "====> Beam is between strips " << iLo << " -- " << iHi << std::endl;
//delete ha;
return;
}
void AnalysisBase::getTDC(){
TProfile *hb = new TProfile("hb","Cluster Charge vs TDC time",12,0,12,100,1000);
Long64_t nentries = fChain->GetEntriesFast();
float lo = 0, hi = 0;
Long64_t nbytes = 0, nb = 0;
std::cout << "============================================" << std::endl;
std::cout << "getTDC(): Determining Optimal TDC time range" << std::endl;
std::cout << "============================================" << std::endl;
for (Long64_t jentry=0; jentry<max(50000,(int)nentries);jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
for(int j=0; j<min((int)clusterNumberPerEvent,10); j++){
//std::cout << "j: " << ", clustersTDC: " << clustersTDC << ", clustersCharge[j]: " << clustersCharge[j] << ", clustersPosition[j]: " << clustersPosition[j] << std::endl;
if(clustersPosition[j] < 0.1) continue;
if(polarity*clustersCharge[j] < kClusterChargeMin) continue;
if(clustersPosition[j]>=iLo && clustersPosition[j]<=iHi && clustersTDC>1.0 &&
polarity*clustersCharge[j]>150 && polarity*clustersCharge[j]<500) hb->Fill(clustersTDC+0.1,polarity*clustersCharge[j]);
}
}
getTDCBins(hb,lo,hi);
tdcLo = lo;
tdcHi = hi;
std::cout << "====> TDC Range determined to be from " << tdcLo << " -- " << tdcHi << std::endl;
hb->Draw();
//delete hb;
return;
}
void AnalysisBase::findBeamRegionAndAlign(int iPass){
cout << "==================================================================" << endl;
cout << "findBeamRegionAndAlign(): Determining Fiducial Region, Pass = " << iPass << endl;
cout << "==================================================================" << endl;
double dxWindow = dxWin;
if(iPass == 1) dxWindow = 1000.0;
if(iPass == 2) dxWindow = 1.0;
TH1F* hthx = new TH1F("hthx","#theta_{X}",1000,-5.0,5.0);
TH1F* hthy = new TH1F("hthy","#theta_{Y}",1000,-5.0,5.0);
TH1F* hx = new TH1F("hx","Y position of matched cluster",800,-40.0,40.0);
TH1F* hxdut = new TH1F("hxdut","Y position of matched cluster",800,-40.0,40.0);
TH1F* hs = new TH1F("hs","Strip number of matched cluster",512,0,512.0);
TH1F* hy = new TH1F("hy","Y position of matched cluster",800,-10.0,10.0);
TH1F* hw = new TH1F("hw","#DeltaX",20000,-100.0,100.0);
TH1F* hn = new TH1F("hn","#DeltaX",4000,-2.0,2.0);
TH1F* hnn = new TH1F("hnn","#DeltaX",400,-0.2,0.2);
TH2F* hxy = new TH2F("hxy","Y_{trk} vs X_{trk}, with cluster",640,-8,8.0,640,-8,8);
TProfile *ht = new TProfile("ht","Cluster Charge vs TDC time",12,0,12,100,1000);
TProfile *hzrot = new TProfile("hzrot","#DeltaX vs Y_{trk} at DUT",1600,-8,8,-1.0,1.0);
TH2F* hca = new TH2F("hca","Y_{trk} vs X_{trk}, with found cluster",160,-8,8.0,320,-8,8);
TH2F* hcf = new TH2F("hcf","Y_{trk} vs X_{trk}, with found cluster",160,-8,8.0,320,-8,8);
hca->Sumw2();
hcf->Sumw2();
double nomStrip=0, detStrip = 0;
Long64_t nbytes = 0, nb = 0;
Int_t nentries = fChain->GetEntriesFast();
for (Long64_t jentry=0; jentry<max(nentries,200000);jentry++) {
Long64_t ientry = LoadTree(jentry);
//if(jentry%1000==0) cout << "At entry = " << jentry << endl;
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
if(n_tp3_tracks > 1) continue;
for(int k=0; k<n_tp3_tracks; k++){
bool goodTime = (clustersTDC >= tdcLo && clustersTDC < tdcHi);
goodTime = clustersTDC>1.0;
if(!goodTime) continue;
double x_trk = vec_trk_tx->at(k)*z_DUT+vec_trk_x->at(k);
double y_trk = vec_trk_ty->at(k)*z_DUT+vec_trk_y->at(k);
transformTrackToDUTFrame(k, x_trk, y_trk, nomStrip, detStrip);
double tx = 1000*vec_trk_tx->at(k);
double ty = 1000*vec_trk_ty->at(k);
double x_trk0 = x_trk;
for(int j=0; j<min(clusterNumberPerEvent,10); j++){
if(clustersPosition[j] < 0.1) continue;
if(polarity*clustersCharge[j] < kClusterChargeMin) continue;
bool goodHit = (clustersPosition[j]>iLo && clustersPosition[j]<iHi);
double x_dut = getDUTHitPosition(j);
x_trk = x_trk0;
double dx = x_dut - x_trk;
hn->Fill(dx);
hnn->Fill(dx);
hca->Fill(x_trk,y_trk);
if(fabs(dx)<dxWindow || iPass==1) {
if(goodHit) {
hx->Fill(x_trk);
hxdut->Fill(x_dut);
hs->Fill(clustersPosition[j]);
hy->Fill(y_trk);
hthx->Fill(tx);
hthy->Fill(ty);
hw->Fill(dx);
hxy->Fill(x_trk,y_trk);
ht->Fill(clustersTDC+0.1,polarity*clustersCharge[j]);
hzrot->Fill(y_trk,dx);
hcf->Fill(x_trk,y_trk);
}
}
}
}
}
//hnn->Draw();
if(iPass==4) {
findCutoutRegion(hcf);
}else {
float xmin=0,xmax=0,ymin=0,ymax=0,txmin=0,txmax=0,tymin=0,tymax=0;
getRange(hthx,txmin,txmax,0.05,2);
getRange(hthy,tymin,tymax,0.05,2);
getRange(hx,xmin,xmax,0.05,2);
getRange(hy,ymin,ymax,0.2,2);
xMin = xmin;
xMax = xmax;
yMin = ymin;
yMax = ymax;
txMin = txmin;
txMax = txmax;
tyMin = tymin;
tyMax = tymax;
float xlo = 0, xhi=0;
if(iPass==1) getRange(hw,xlo,xhi,0.1,1);
if(iPass==2) getRange(hn,xlo,xhi,0.1,1);
if(iPass==3) getRange(hnn,xlo,xhi,0.1,1);
double XOFF = (xhi + xlo)/2.0;
xOff = xOff - XOFF;
double x_ave = (xMax + xMin)/2.0;
double y_ave = (yMax + yMin)/2.0;
xGloOff = xGloOff - x_ave;
yGloOff = yGloOff - y_ave;
// Check/correct for z rotation
if(iPass==3 && correctForZRotation){
cout << "Checking for z-rotation" << endl;
TF1* p1 = new TF1("p1","[0]+[1]*x",yMin,yMax);
p1->SetParameters(0.0,0.0001);
hzrot->Fit(p1,"0R");
double rz = p1->GetParameter(1);
cout << "=======================================================================" << endl;
cout << "==> Updating z rotation angle from " << Rz << " to " << Rz-rz << " mrad" << endl;
Rz = Rz - rz;
delete p1;
}
float lo = 0, hi = 0;
getTDCBins(ht,lo,hi);
tdcLo = lo;
tdcHi = hi;
cout << "=================================================================" << endl;
cout << "====> TDC Range updated to be from " << tdcLo << " -- " << tdcHi << std::endl;
cout << "=================================================================" << endl;
cout << "====> Fiducial regions, xLo, xHi (Strip numbers) = " << iLo << " " << iHi << endl;
cout << "====> Fiducial regions, xLo, xHi (pos in mm) = " << xMin << " " << xMax << " mm " << endl;
cout << "====> Fiducial regions, yLo, yHi (pos in mm) = " << yMin << " " << yMax << " mm " << endl;
cout << "====> Fiducial regions, thxLo, thxHi (pos in mrad) = " << txMin << " " << txMax << " mrad" << endl;
cout << "====> Fiducial regions, thyLo, thxHi (pos in mrad) = " << tyMin << " " << tyMax << " mrad" << endl;
cout << "=================================================================" << endl;
cout << "====> Shifting sensor by dX = " << XOFF << " mm " << ", new xOff = " << xOff << endl;
cout << "====> Global X, Y shifts: " << x_ave << " " << y_ave << endl;
cout << "====> New global x,y = " << xGloOff << " " << yGloOff << endl;
}
//if(iPass==1) hw->Draw(); // for debugging
//if(iPass==2) hy->Draw();
//return;
delete hxdut;
delete hthx;
delete hthy;
delete hx;
delete hy;
delete hxy;
delete hw;
delete hn;
delete hnn;
delete hs;
delete ht;
delete hzrot;
delete hcf;
delete hca;
return;
}
Double_t AnalysisBase::getCorrChannel(double ch){
//if(ch>=128 && ch<=255) return (ch + stripGap[0]);
//if(ch>=256 && ch<=383) return (ch + stripGap[1]);
//if(ch>=384 && ch<=512) return (ch + stripGap[2]);
if(ch<=128) return (ch - stripGap[0]);
if(ch>=128 && ch<=255) return (ch - stripGap[1]);
if(ch>=256 && ch<=383) return (ch - stripGap[2]);
if(ch>=384 && ch<=512) return (ch - stripGap[3]);
return ch;
}
void AnalysisBase::findChipBoundary(){
// Look for holes in hit profile, due to boundary between chips
std::cout << "=======================================================================================" << std::endl;
std::cout << "findChipBoundary(): Looking for gaps due to unconnected strips in between Beetle chips " << std::endl;
std::cout << "=======================================================================================" << std::endl;
TH1F* hx2 = new TH1F("hx2","X position of cluster",200,-10.0,10.0);
Long64_t nbytes = 0, nb = 0;
Long64_t nentries = fChain->GetEntriesFast();
for (Long64_t jentry=0; jentry<max((int)nentries,20000);jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
for(int j=0; j<min((int)clusterNumberPerEvent,10); j++){
if(clustersPosition[j] < 0.1) continue;
if(polarity*clustersCharge[j] < kClusterChargeMin) continue;
bool goodHit = (clustersPosition[j]>iLo && clustersPosition[j]<iHi);
goodHit = goodHit && (clustersTDC >= tdcLo && clustersTDC <= tdcHi);
double x_dut = getDUTHitPosition(j);
hx2->Fill(x_dut);
}
}
//hx2->Draw(); // for debugging only
int i1 = hx2->FindBin(xMin);
int i2 = hx2->FindBin(xMax);
double ave = 0.0;
int ixL=i1, ixH=i2;
for(int i=i1; i<i2; i++){
ave = ave + hx2->GetBinContent(i);
}
ave = ave / (i2 - i1 + 1);
for(int i=i1; i<i2; i++){
double v = hx2->GetBinContent(i);
double xx = hx2->GetBinCenter(i);
double r = v/ave;
if(r < ratioForHole && xx < xLeftHole) {
xLeftHole = xx;
ixL = i;
}
if(r < ratioForHole && xx > xRightHole) {
xRightHole = xx;
ixH = i;
}
}
// Count the number of "dead" channels in this region
float ndead = 0;
for(int i=ixL; i<=ixH; i++){
double v = hx2->GetBinContent(i);
double r = v/ave;
if(r < ratioForHole) ndead = ndead + 1;
}
double frdead = ndead / (ixH - ixL + 1);
if(frdead>0.8 && ndead>4){
xLeftHole = xLeftHole - 0.2;
xRightHole = xRightHole + 0.2;
}else{
xLeftHole = 999;
xRightHole = -999;
}
hx2->Draw();
//std::cout << "Xmin, xMax = " << xMin << " " << xMax << std::endl;
//std::cout << "Bin Edges: " << i1 << " " << i2 << std::endl;
//std::cout << "Hole position: " << xLeftHole << " " << xRightHole << " , fraction dead = " << ndead << " " << frdead << std::endl;
delete hx2;
}
void AnalysisBase::findDeadRegions(){
cout << "=================================================" << endl;
cout << "findDeadRegions(): Looking for dead strip regions " << endl;
cout << "=================================================" << endl;
TH1F* hf = new TH1F("hf","X position of matched cluster",200,-10.0,10.0);
TH1F* hnf = new TH1F("hnf","X position of matched cluster",200,-10.0,10.0);
double nomStrip = 0, detStrip = 0;
Long64_t nbytes = 0, nb = 0;
Int_t nentries = fChain->GetEntriesFast();
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
if(n_tp3_tracks > 1) continue;
// Loop over TPIX tracks in event
for(int k=0; k<n_tp3_tracks; k++){
double x_trk = vec_trk_tx->at(k)*z_DUT+vec_trk_x->at(k);
double y_trk = vec_trk_ty->at(k)*z_DUT+vec_trk_y->at(k);
transformTrackToDUTFrame(k, x_trk, y_trk, nomStrip, detStrip );
bool goodTrack = false;
bool inFiducial = false;
if(x_trk>xMin && x_trk<xMax && y_trk>yMin && y_trk<yMax) inFiducial = true;
inFiducial = inFiducial && (x_trk<xLeftHole || x_trk>xRightHole);
double tx = 1000*vec_trk_tx->at(k);
double ty = 1000*vec_trk_ty->at(k);
if(tx>txMin && tx<txMax && ty>tyMin && ty<tyMax) goodTrack = true;
bool goodTime = (clustersTDC >= tdcLo && clustersTDC <= tdcHi);
bool foundHit = false;
// Loop over clusters
for(int j=0; j<min(clusterNumberPerEvent,10); j++){
if(clustersPosition[j] < 0.1) continue;
if(polarity*clustersCharge[j] < kClusterChargeMin) continue;
//bool goodHit = (clustersPosition[j]>iLo || clustersPosition[j]<iHi);
double x_dut = getDUTHitPosition(j);
double dx = x_dut - x_trk;
if(goodTrack && inFiducial && goodTime && fabs(dx)<dxWin) {
foundHit = true;
}
}
if(inFiducial && goodTrack && goodTime) {
hf->Fill(x_trk);
if(!foundHit) hnf->Fill(x_trk);
}
}
}
TH1F *hineff = (TH1F*)hnf->Clone("hineff");
hineff->Divide(hf);
// Only for debugging
/*
TCanvas *cc = new TCanvas("cc","Hits",800,800);
cc->Divide(2,2);
cc->cd(1);
hnf->Draw();
cc->cd(2);
hf->Draw();
cc->cd(3);
hineff->Draw();
*/
for(int k = 0; k<hineff->GetNbinsX()-1; k++){
if(nDeadRegion > 20){
cout << "WARNING: Exceeded 20 dead regions, exiting from loop" << endl;
break;
}
double v1 = hineff->GetBinContent(k);
if(v1 > k_DeadStrip){
deadRegionLo[nDeadRegion] = hf->GetBinLowEdge(k);
for(int j = k; j<hineff->GetNbinsX(); j++){
double v2 = hineff->GetBinContent(j);
if(v2 < k_DeadStrip){
deadRegionHi[nDeadRegion] = hf->GetBinLowEdge(j);
cout << "====> Found inefficiency strip region from x : "
<< deadRegionLo[nDeadRegion] << " <===> " << deadRegionHi[nDeadRegion] << " mm" << endl;
k = j + 1;
nDeadRegion++;
break;
}
}
}
}
delete hf;
delete hnf;
delete hineff;
return;
}
void AnalysisBase::transformTrackToDUTFrame(int k, double& x_trk, double& y_trk, double& elecStrip, double& detStrip){
y_trk = y_trk + yGloOff;
double dzs = x_trk * sin(Ry);
x_trk = x_trk - Rz*y_trk;
y_trk = y_trk + Rz*x_trk;
x_trk = x_trk + vec_trk_tx->at(k)*dzs;
x_trk = x_trk/cos(Ry);
double rStrip = nChan - (x_trk-xOff)/stripPitch;
double corStrip = getCorrChannel(rStrip);
//double dx = (corStrip - rStrip)*stripPitch;
//x_trk = x_trk + dx;
detStrip = rStrip; //nChan - (x_trk-xOff)/stripPitch;
elecStrip = corStrip;
x_trk = x_trk + xGloOff;
return;
}
void AnalysisBase::PrepareDUT(){
//----------------------------
// Get Beam Location (Strips)
//----------------------------
getBeamLoc();
//return;
// ==================================
// Determining Optimal TDC time range
// ==================================
getTDC();
// --------------------------------------------------------------
// Find beam fiducial region (slopes & Y range) & align residuals
// ---------------------------------------------------------------
findBeamRegionAndAlign(1);
findBeamRegionAndAlign(2);
yMid = yMin + (yMax-yMin)/2.0;
yHi2 = yMax - 2.0;
findBeamRegionAndAlign(3);
yMid = yMin + (yMax-yMin)/2.0;
yHi2 = yMax - 2.0;
//return;
if(yInt2[1] > yMax){
yInt1[0] = yMax-0.8; yInt1[1] = yMax-0.4;
yInt2[0] = yMax-0.4; yInt2[1] = yMax;
yInt3[0] = yMin; yInt3[1] = yMax-0.8;
}
cout << "Y Region Definitions:" << endl;
cout << " ===> Region 1: " << yInt1[0] << " -- " << yInt1[1] << " mm" << endl;
cout << " ===> Region 2: " << yInt2[0] << " -- " << yInt2[1] << " mm" << endl;
cout << " ===> Region 3: " << yInt3[0] << " -- " << yInt3[1] << " mm" << endl;
//return;
// Correct for gaps between Beetle chips
correctForStripGaps();
//findBeamRegionAndAlign(4);
if(holeSector) findBeamRegionAndAlign(4);
findChipBoundary();
cout << "====> Hole position: " << xLeftHole << " " << xRightHole << endl;
//return;
if(vetoDeadStripRegions) findDeadRegions();
//return;
setCrossTalkCorr();
}
Double_t AnalysisBase::getDUTHitPosition(int j){
double pos = getCorrChannel(clustersPosition[j]);
//double pos = clustersPosition[j];
double x_dut = (nChan - pos)*stripPitch + xOff;
x_dut = x_dut + xGloOff;
return x_dut;
}
void AnalysisBase::setCrossTalkCorr(){
float biasVal = atof(m_bias);
chargeCorrSlopeOdd = 0.0;
chargeCorrSlopeEven = 0.0;
cout << "Bias Value = " << m_bias << " " << biasVal << endl;
if(m_board.Contains("A2") && m_sector=="1" && biasVal < 260){
chargeCorrSlopeOdd = 0.435;
chargeCorrSlopeEven = 0.400;
}else if(m_board.Contains("A2") && m_sector=="1" && biasVal > 320){
chargeCorrSlopeOdd = 0.148;
chargeCorrSlopeEven = 0.094;
}else if(m_board.Contains("A2") && m_sector=="1" && biasVal == 300){
chargeCorrSlopeOdd = 0.130;
chargeCorrSlopeEven = 0.088;
}else if(m_board.Contains("A2") && m_sector=="2"){
chargeCorrSlopeOdd = 0.120;
chargeCorrSlopeEven = 0.085;
}
}
void AnalysisBase::findCutoutRegion(TH2F* h){
cout << "===============================================================" << endl;
cout << "findCutoutRegion(): Looking for cutout region in D type sensor " << endl;
cout << "===============================================================" << endl;
int ilx = h->GetXaxis()->FindBin(xMin) + 1;
int ihx = h->GetXaxis()->FindBin(xMax) - 1;
int ily = 1;//h->GetYaxis()->GetBinFindBin(yMin) - 10;
int ihy = h->GetYaxis()->FindBin(yMax) + 10;
if(ily < 0) ily = 0;
if(ihy > h->GetYaxis()->GetNbins()) ihy = h->GetYaxis()->GetNbins();
int nb = h->GetXaxis()->GetNbins();
double xlow = h->GetXaxis()->GetBinLowEdge(1);
double xhi = h->GetXaxis()->GetBinLowEdge(nb)+h->GetXaxis()->GetBinWidth(1);
TH1F *hpr = new TH1F("hpr","Y vs X, edge",nb,xlow,xhi);
TH1D *hpy;
int ipeak = 0;
for(int i=ilx; i<ihx-1; i++){
hpy = h->ProjectionY("hpy",i,i);
double maxcon = 0.0;
double maxbin = 0.0;
for(int j=ily;j<ihy;j++){
maxcon = maxcon + hpy->GetBinContent(j);
if(hpy->GetBinContent(j) > maxbin) {
maxbin = hpy->GetBinContent(j) ;
ipeak = j;
}
}
// Find the lower "edge"
for(int j=ipeak; j>=ily; j--){
double r0 = (hpy->GetBinContent(j-1)+hpy->GetBinContent(j-2))/2.0;
double r1 = hpy->GetBinContent(j);
double r2 = (hpy->GetBinContent(j+3)+hpy->GetBinContent(j+4))/2.0;
if(r1<=2 && r0<=1.5 && r2>5*r1 && r2>8){
hpr->SetBinContent(i,hpy->GetBinCenter(j)+2*hpy->GetBinWidth(j));
hpr->SetBinError(i,hpy->GetBinWidth(1)/2.0);
//cout << "found: " << i << " " << hpy->GetEntries() << " " << hpy->GetBinCenter(j) << " "
// << hpy->GetBinContent(j) << " " << r0 << " " << r1 << " " << r2 << " " << endl;
break;
}
}
}
TF1* poly2 = new TF1("poly2","[0]+[1]*x+[2]*x*x",xMin,xMax);
poly2->SetParameters(-1.5,-0.17,-0.15);
hpr->Fit(poly2,"R0");
hpr->SetLineColor(kRed);
holeQuadPar[0] = poly2->GetParameter(0);
holeQuadPar[1] = poly2->GetParameter(1);
holeQuadPar[2] = poly2->GetParameter(2);
//return;
delete poly2;
delete hpr;
delete hpy;
return;
}
bool AnalysisBase::isInCutoutRegion(double xtrk, double ytrk){
// Some protections here..
if(!removeTracksInHole) return false;
if(holeQuadPar[0]==0 || holeQuadPar[1]==0 || holeQuadPar[2]==0) return false;
// Ok, looks like we mean to really remove these tracks
double yhole = holeQuadPar[0]+holeQuadPar[1]*xtrk+holeQuadPar[2]*xtrk*xtrk;
if(ytrk < yhole + minDistFromHole) return true;
return false;
}
double AnalysisBase::DistToCutoutRegion(double xtrk, double ytrk){
// Some protections here..
if(holeQuadPar[0]==0 || holeQuadPar[1]==0 || holeQuadPar[2]==0) return 999.0;
// Ok, looks like we mean to really remove these tracks
double yhole = holeQuadPar[0]+holeQuadPar[1]*xtrk+holeQuadPar[2]*xtrk*xtrk;
double a = 2.0*holeQuadPar[2]*xtrk+holeQuadPar[1];
double c = yhole - a*xtrk;
double b = -1.0;
double dist = fabs(a*xtrk + b*ytrk + c) / sqrt(a*a+b*b);
if(ytrk < yhole) dist = -1.0*dist;
return dist;
}
void AnalysisBase::correctForStripGaps(){
cout << "================================================================================" << endl;
cout << "correctForStripGaps(): Correcting x fiducial range for gaps between Beetle chips" << endl;
cout << "================================================================================" << endl;
double x1 = getCorrChannel(iLo);
double x2 = getCorrChannel(iHi);
xMax = (nChan - x1)*stripPitch + xOff + xGloOff;
xMin = (nChan - x2)*stripPitch + xOff + xGloOff;
double xave = (xMax + xMin)/2.;
xGloOff = xGloOff - xave;
xMin = xMin - xave;
xMax = xMax - xave;
std::cout << "====> Updating Global Offset to use Strips, *** New xMin, xMax == " << xMin << " " << xMax << endl;
std::cout << "====> xMin, xMax = " << xMin << " " << xMax << endl;
std::cout << "====> xMin, xMax = " << yMin << " " << yMax << endl;
return;
}
iaro