Newer
Older
TB_Chris / TbUT / scripts / .svn / text-base / AnalysisBase.C.svn-base
  1. //#define AnalysisBase_cxx
  2. #include "AnalysisBase.h"
  3.  
  4.  
  5. AnalysisBase::AnalysisBase(TTree *tree) : fChain(0)
  6. {
  7.  
  8. // Notes
  9. // For Rz alignment, one should SUBTRACT the slope found from the deltaX vs Y plot.
  10. // For Z alignment, one should ADD the slope found from the deltaX vs Y plot.
  11. // For Ry: If slope of DX vs X is negative (positive), reduce (increase) Ry
  12. //
  13. //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
  14.  
  15. minDistFromHole = m_minDistFromHole;
  16. holeSector = m_board.Contains("D") && (m_sector=="1" || m_sector=="2" || m_sector == "3");
  17. if(!holeSector) minDistFromHole = -9999;
  18. stripPitch = 0.190;
  19. z_DUT = 370; //This needs to be set for different boards
  20. Rz = -0.0218;
  21. Ry = 0.00;
  22. dxWin = 0.25;
  23. xGloOff = -9.2;
  24. yGloOff = -7.5;
  25. xOff = -54.78;
  26. correctForZRotation = true;
  27. channelOffset = 0.0;
  28. xLeftHole = 999;
  29. xRightHole = -999;
  30. stripGap[0] = 0;
  31. stripGap[1] = 0;
  32. stripGap[2] = 0;
  33. stripGap[3] = 0;
  34. nDeadRegion = 0;
  35. if( m_board.Contains("D")) {
  36. stripPitch = 0.095;
  37. channelOffset = 512.0;
  38. //stripGap[0]=19;
  39. //stripGap[1]=39;
  40. //stripGap[2]=58;
  41. //stripGap[3]=0;
  42.  
  43. stripGap[0]=58;
  44. stripGap[1]=39;
  45. stripGap[2]=19;
  46. stripGap[3]=0;
  47.  
  48. if(m_board.Contains("D5")) {
  49. z_DUT = 1100.0;
  50. //Rz = -0.0230;
  51. }else if(m_board.Contains("D7")) {
  52. z_DUT = 1123.0;
  53. if(m_sector=="5") z_DUT = 1087.0;
  54. if(m_sector=="1") z_DUT = 1040.0;
  55. //Ry = -0.008;
  56. }
  57. }
  58.  
  59.  
  60.  
  61. yInt1[0] = 2.8; yInt1[1] = 3.4;
  62. yInt2[0] = 3.4; yInt2[1] = 4.1;
  63. yInt3[0] = 1.0; yInt3[1] = 2.3;
  64.  
  65. if(m_board.Contains("A8")){
  66. z_DUT = 1100.0;
  67. Rz = -0.0286;
  68. yInt1[0] = 1.8; yInt1[1] = 2.0;
  69. yInt2[0] = 2.0; yInt2[1] = 2.4;
  70. yInt3[0] = 1.0; yInt3[1] = 2.3;
  71. if(m_sector=="3" || m_sector=="6") {
  72. channelOffset = -128.0*2+22;
  73. }else{
  74. channelOffset = -128.0;
  75. }
  76. if(m_scanType.Contains("Angle")){
  77. if(m_angle == "10") {
  78. Ry = 9.3*(TMath::Pi()/180.0);
  79. z_DUT = z_DUT - 18.0;
  80. } else if(m_angle == "20") {
  81. Ry = 19*(TMath::Pi()/180.0);
  82. z_DUT = z_DUT - 18.0;
  83. }
  84. }
  85. }else if( m_board.Contains("A4")){
  86. z_DUT = 381.0;
  87. Rz = -0.0238;
  88. if(m_sector=="3" || m_sector=="6") {
  89. channelOffset = -128.0*2.0+19;
  90. }else{
  91. channelOffset = -128.0;
  92. }
  93. }else if( m_board.Contains("A6")){
  94. yInt1[0] = 2.5; yInt1[1] = 3.1;
  95. yInt2[0] = 3.1; yInt2[1] = 3.7;
  96. yInt3[0] = 1.0; yInt3[1] = 2.3;
  97. if(m_sector=="3" || m_sector=="6") {
  98. channelOffset = -128.0*2+22;
  99. }else{
  100. channelOffset = -128.0;
  101. }
  102. }else if( m_board.Contains("A2") || m_board.Contains("A1")){
  103. z_DUT = 310;
  104. yInt1[0] = -5.0; yInt1[1] = -2.0;
  105. yInt2[0] = -2.0; yInt2[1] = 1.0;
  106. yInt3[0] = 1.0; yInt3[1] = 5.0;
  107. }
  108.  
  109. // Board A1 alignment is really messed up for a few runs! ADHOC here!
  110. if(m_board.Contains("A1")){// && (m_bias=="300" || m_bias=="350" || m_bias=="150" || m_bias==")){
  111. dxWin = 1.0;
  112. correctForZRotation = false;
  113. }
  114. holeQuadPar[0] = 0;
  115. holeQuadPar[1] = 0;
  116. holeQuadPar[2] = 0;
  117. removeTracksInHole = removeTracksInHoleDef;
  118. // Only look to remove tacks from hole area if D-type and in sectors 1, 2, or 3.
  119. if(!m_board.Contains("D") || !(m_sector=="1" || m_sector=="2" || m_sector=="3") ) removeTracksInHole = false;
  120.  
  121. polarity = 1.0;
  122. if(isPType) polarity = -1;
  123.  
  124.  
  125. }
  126.  
  127. AnalysisBase::~AnalysisBase()
  128. {
  129. if (!fChain) return;
  130. delete fChain->GetCurrentFile();
  131. }
  132.  
  133. Int_t AnalysisBase::GetEntry(Long64_t entry)
  134. {
  135. // Read contents of entry.
  136. if (!fChain) return 0;
  137. return fChain->GetEntry(entry);
  138. }
  139. Long64_t AnalysisBase::LoadTree(Long64_t entry)
  140. {
  141. // Set the environment to read one entry
  142. if (!fChain) return -5;
  143. Long64_t centry = fChain->LoadTree(entry);
  144. if (centry < 0) return centry;
  145. if (fChain->GetTreeNumber() != fCurrent) {
  146. fCurrent = fChain->GetTreeNumber();
  147. Notify();
  148. }
  149. return centry;
  150. }
  151.  
  152. void AnalysisBase::Init(TTree *tree)
  153. {
  154. // The Init() function is called when the selector needs to initialize
  155. // a new tree or chain. Typically here the branch addresses and branch
  156. // pointers of the tree will be set.
  157. // It is normally not necessary to make changes to the generated
  158. // code, but the routine can be extended by the user if needed.
  159. // Init() will be called many times when running on PROOF
  160. // (once per file to be processed).
  161.  
  162. // Set object pointer
  163. vec_trk_x = 0;
  164. vec_trk_y = 0;
  165. vec_trk_tx = 0;
  166. vec_trk_ty = 0;
  167. vec_trk_chi2ndf = 0;
  168. // Set branch addresses and branch pointers
  169. if (!tree) return;
  170. fChain = tree;
  171. fCurrent = -1;
  172. fChain->SetMakeClass(1);
  173.  
  174. fChain->SetBranchAddress("clusterNumberPerEvent", &clusterNumberPerEvent, &b_clusterNumberPerEvent);
  175. fChain->SetBranchAddress("clustersTDC", &clustersTDC, &b_clustersTDC);
  176. fChain->SetBranchAddress("timestamps", &timestamps, &b_timestamps);
  177. fChain->SetBranchAddress("clustersPosition", clustersPosition, &b_clustersPosition);
  178. fChain->SetBranchAddress("clustersSeedPosition", clustersSeedPosition, &b_clustersSeedPosition);
  179. fChain->SetBranchAddress("clustersCharge", clustersCharge, &b_clustersCharge);
  180. fChain->SetBranchAddress("clustersSize", clustersSize, &b_clustersSize);
  181. fChain->SetBranchAddress("clustersSeedCharge", clustersSeedCharge, &b_clustersSeedCharge);
  182. fChain->SetBranchAddress("clustersCharge2StripLeft", clustersCharge2StripLeft, &b_clustersCharge2StripLeft);
  183. fChain->SetBranchAddress("clustersCharge1StripLeft", clustersCharge1StripLeft, &b_clustersCharge1StripLeft);
  184. fChain->SetBranchAddress("clustersCharge1StripRight", clustersCharge1StripRight, &b_clustersCharge1StripRight);
  185. fChain->SetBranchAddress("clustersCharge2StripRight", clustersCharge2StripRight, &b_clustersCharge2StripRight);
  186. fChain->SetBranchAddress("n_tp3_tracks", &n_tp3_tracks, &b_n_tp3_tracks);
  187. fChain->SetBranchAddress("vec_trk_x", &vec_trk_x, &b_vec_trk_x);
  188. fChain->SetBranchAddress("vec_trk_y", &vec_trk_y, &b_vec_trk_y);
  189. fChain->SetBranchAddress("vec_trk_tx", &vec_trk_tx, &b_vec_trk_tx);
  190. fChain->SetBranchAddress("vec_trk_ty", &vec_trk_ty, &b_vec_trk_ty);
  191. fChain->SetBranchAddress("vec_trk_chi2ndf", &vec_trk_chi2ndf, &b_vec_trk_chi2ndf);
  192. fChain->SetBranchAddress("dtime", &dtime, &b_dtime);
  193.  
  194. Notify();
  195. }
  196.  
  197. Bool_t AnalysisBase::Notify()
  198. {
  199. // The Notify() function is called when a new file is opened. This
  200. // can be either for a new TTree in a TChain or when when a new TTree
  201. // is started when using PROOF. It is normally not necessary to make changes
  202. // to the generated code, but the routine can be extended by the
  203. // user if needed. The return value is currently not used.
  204.  
  205. return kTRUE;
  206. }
  207.  
  208. void AnalysisBase::Show(Long64_t entry)
  209. {
  210. // Print contents of entry.
  211. // If entry is not specified, print current entry
  212. if (!fChain) return;
  213. fChain->Show(entry);
  214. }
  215. /*
  216. Int_t AnalysisBase::Cut(Long64_t entry)
  217. {
  218. // This function may be called from Loop.
  219. // returns 1 if entry is accepted.
  220. // returns -1 otherwise.
  221. return 1;
  222. }
  223. */
  224. void AnalysisBase::getRange(TH1 *h, float &lo, float& hi, float thresh, int nSkipMax){
  225.  
  226. int nfail = 0;
  227. int iMaxBin = h->GetMaximumBin();
  228. lo = iMaxBin;
  229. hi = iMaxBin;
  230. double maxValue = h->GetBinContent(iMaxBin);
  231. // Only for debugging....
  232. //h->Draw();
  233. //std::cout << "maxValue = " << iMaxBin << " " << maxValue << std::endl;
  234.  
  235. for(int j=iMaxBin;j>=0;j--){
  236. double v = h->GetBinContent(j);
  237. double r = v / maxValue;
  238. //std::cout << "r = " << r << std::endl; // Only for debugging
  239. if( r > thresh) {
  240. lo = j ;
  241. nfail = 0;
  242. }else if( r <= thresh) {
  243. nfail++;
  244. }
  245. if(nfail >= nSkipMax) break;
  246. }
  247.  
  248. nfail = 0;
  249. for(int j=iMaxBin;j<=h->GetNbinsX();j++){
  250. double v = h->GetBinContent(j);
  251. double r = v / maxValue;
  252. //std::cout << "r = " << r << std::endl; // Only for debugging
  253. if( r > thresh) {
  254. hi = j;
  255. nfail = 0;
  256. }else if( r <= thresh) {
  257. nfail++;
  258. }
  259. if(nfail >= nSkipMax) break;
  260. }
  261. lo = h->GetBinCenter(lo);
  262. hi = h->GetBinCenter(hi);
  263.  
  264. //std::cout << "lo, hi = " << lo << " " << hi << std::endl; // only for debugging
  265.  
  266. return;
  267. }
  268.  
  269.  
  270. void AnalysisBase::getTDCBins(TProfile* h, float& lo, float& hi){
  271. getRange(h,lo,hi,0.97,1);
  272. lo = lo - 1.0 + 0.5;
  273. hi = hi + 0.5;
  274. return;
  275. }
  276.  
  277.  
  278.  
  279. void AnalysisBase::getBeamLocation(TH1F *h, float &lo, float& hi){
  280. //double dif = 10;
  281. for(int j = 0; j<10; j++){
  282. std::cout << "====> Finding Beam, iteration " << j+1 << std::endl;
  283. getRange(h,lo,hi);
  284. int dif = hi - lo + 1;
  285. if(dif < 5){
  286. int i1 = h->FindBin(lo);
  287. int i2 = h->FindBin(hi);
  288. for(int k=i1; k<=i2; k++){
  289. h->SetBinContent(k,0);
  290. }
  291. }else{
  292. break;
  293. }
  294. }
  295.  
  296. //h->Draw();
  297. return;
  298. }
  299.  
  300. double AnalysisBase::getXOffset(TH1F *h1w){
  301. float xlo = 0, xhi=0;
  302. getRange(h1w,xlo,xhi,0.1,1);
  303. return (xlo+xhi)/2.0;
  304.  
  305. }
  306.  
  307. void AnalysisBase::getBeamLoc(){
  308.  
  309. TH1F* ha = new TH1F("ha","Strip # of cluster with track",512,0.0,512);
  310. //TH1F* hb = new TH1F("hb","Strip # of cluster with track",512,0.0,512);
  311. Long64_t nentries = fChain->GetEntriesFast();
  312. float lo = 0, hi = 0;
  313.  
  314. Long64_t nbytes = 0, nb = 0;
  315. std::cout << "======================================= " << std::endl;
  316. std::cout << "getBeamLoc(): Determining Beam Position " << nentries << std::endl;
  317. std::cout << "======================================= " << std::endl;
  318. for (Long64_t jentry=0; jentry<max(50000,(int)nentries);jentry++) {
  319. Long64_t ientry = LoadTree(jentry);
  320. if (ientry < 0) break;
  321. nb = fChain->GetEntry(jentry); nbytes += nb;
  322. //cout << clusterNumberPerEvent << endl;
  323. for(int j=0; j<min((int)clusterNumberPerEvent,10); j++){
  324. //std::cout << clustersPosition[j] << std::endl;
  325. //cout << polarity << " " << polarity*clustersCharge[j] << endl;
  326. if(polarity*clustersCharge[j] < kClusterChargeMin) continue;
  327. int iChan = clustersSeedPosition[j];
  328. //if(j<500) cout << "iChan = " << iChan << clustersCharge[j] << " " << 5*noise[iChan] << endl;
  329. if(polarity*clustersCharge[j]<4*noise[iChan]) continue;
  330. if(clustersPosition[j]>0.1&&clustersSize[j]==1) ha->Fill(clustersPosition[j]);
  331. if(clustersPosition[j]>0.1&&clustersSize[j]==2) ha->Fill(clustersPosition[j]);
  332. }
  333. }
  334. int num = ha->GetEntries();
  335. for(int i=0; i<1000; i++){
  336. if(num < 1000) continue;
  337. break;
  338. }
  339.  
  340. if(num < 1000){
  341. std::cout << "ERROR: Something wrong here, insufficient entries in GetBeamLoc(), nEntries = " << num << std::endl;
  342. //exit(1);
  343. }
  344. //ha->Draw();
  345. getBeamLocation(ha,lo,hi);
  346. iLo = lo;
  347. iHi = hi;
  348. std::cout << "====> Beam is between strips " << iLo << " -- " << iHi << std::endl;
  349.  
  350. //delete ha;
  351. return;
  352.  
  353. }
  354.  
  355. void AnalysisBase::getTDC(){
  356. TProfile *hb = new TProfile("hb","Cluster Charge vs TDC time",12,0,12,100,1000);
  357. Long64_t nentries = fChain->GetEntriesFast();
  358. float lo = 0, hi = 0;
  359.  
  360. Long64_t nbytes = 0, nb = 0;
  361. std::cout << "============================================" << std::endl;
  362. std::cout << "getTDC(): Determining Optimal TDC time range" << std::endl;
  363. std::cout << "============================================" << std::endl;
  364. for (Long64_t jentry=0; jentry<max(50000,(int)nentries);jentry++) {
  365. Long64_t ientry = LoadTree(jentry);
  366. if (ientry < 0) break;
  367. nb = fChain->GetEntry(jentry); nbytes += nb;
  368.  
  369. for(int j=0; j<min((int)clusterNumberPerEvent,10); j++){
  370. //std::cout << "j: " << ", clustersTDC: " << clustersTDC << ", clustersCharge[j]: " << clustersCharge[j] << ", clustersPosition[j]: " << clustersPosition[j] << std::endl;
  371. if(clustersPosition[j] < 0.1) continue;
  372. if(polarity*clustersCharge[j] < kClusterChargeMin) continue;
  373. if(clustersPosition[j]>=iLo && clustersPosition[j]<=iHi && clustersTDC>1.0 &&
  374. polarity*clustersCharge[j]>150 && polarity*clustersCharge[j]<500) hb->Fill(clustersTDC+0.1,polarity*clustersCharge[j]);
  375. }
  376. }
  377.  
  378. getTDCBins(hb,lo,hi);
  379. tdcLo = lo;
  380. tdcHi = hi;
  381. std::cout << "====> TDC Range determined to be from " << tdcLo << " -- " << tdcHi << std::endl;
  382. hb->Draw();
  383. //delete hb;
  384.  
  385. return;
  386. }
  387.  
  388. void AnalysisBase::findBeamRegionAndAlign(int iPass){
  389.  
  390. cout << "==================================================================" << endl;
  391. cout << "findBeamRegionAndAlign(): Determining Fiducial Region, Pass = " << iPass << endl;
  392. cout << "==================================================================" << endl;
  393.  
  394. double dxWindow = dxWin;
  395. if(iPass == 1) dxWindow = 1000.0;
  396. if(iPass == 2) dxWindow = 1.0;
  397.  
  398. TH1F* hthx = new TH1F("hthx","#theta_{X}",1000,-5.0,5.0);
  399. TH1F* hthy = new TH1F("hthy","#theta_{Y}",1000,-5.0,5.0);
  400. TH1F* hx = new TH1F("hx","Y position of matched cluster",800,-40.0,40.0);
  401. TH1F* hxdut = new TH1F("hxdut","Y position of matched cluster",800,-40.0,40.0);
  402. TH1F* hs = new TH1F("hs","Strip number of matched cluster",512,0,512.0);
  403. TH1F* hy = new TH1F("hy","Y position of matched cluster",800,-10.0,10.0);
  404. TH1F* hw = new TH1F("hw","#DeltaX",20000,-100.0,100.0);
  405. TH1F* hn = new TH1F("hn","#DeltaX",4000,-2.0,2.0);
  406. TH1F* hnn = new TH1F("hnn","#DeltaX",400,-0.2,0.2);
  407. TH2F* hxy = new TH2F("hxy","Y_{trk} vs X_{trk}, with cluster",640,-8,8.0,640,-8,8);
  408. TProfile *ht = new TProfile("ht","Cluster Charge vs TDC time",12,0,12,100,1000);
  409. TProfile *hzrot = new TProfile("hzrot","#DeltaX vs Y_{trk} at DUT",1600,-8,8,-1.0,1.0);
  410. TH2F* hca = new TH2F("hca","Y_{trk} vs X_{trk}, with found cluster",160,-8,8.0,320,-8,8);
  411. TH2F* hcf = new TH2F("hcf","Y_{trk} vs X_{trk}, with found cluster",160,-8,8.0,320,-8,8);
  412. hca->Sumw2();
  413. hcf->Sumw2();
  414.  
  415. double nomStrip=0, detStrip = 0;
  416. Long64_t nbytes = 0, nb = 0;
  417. Int_t nentries = fChain->GetEntriesFast();
  418. for (Long64_t jentry=0; jentry<max(nentries,200000);jentry++) {
  419. Long64_t ientry = LoadTree(jentry);
  420. //if(jentry%1000==0) cout << "At entry = " << jentry << endl;
  421. if (ientry < 0) break;
  422. nb = fChain->GetEntry(jentry); nbytes += nb;
  423. if(n_tp3_tracks > 1) continue;
  424.  
  425. for(int k=0; k<n_tp3_tracks; k++){
  426. bool goodTime = (clustersTDC >= tdcLo && clustersTDC < tdcHi);
  427. goodTime = clustersTDC>1.0;
  428. if(!goodTime) continue;
  429. double x_trk = vec_trk_tx->at(k)*z_DUT+vec_trk_x->at(k);
  430. double y_trk = vec_trk_ty->at(k)*z_DUT+vec_trk_y->at(k);
  431.  
  432. transformTrackToDUTFrame(k, x_trk, y_trk, nomStrip, detStrip);
  433.  
  434. double tx = 1000*vec_trk_tx->at(k);
  435. double ty = 1000*vec_trk_ty->at(k);
  436.  
  437. double x_trk0 = x_trk;
  438.  
  439. for(int j=0; j<min(clusterNumberPerEvent,10); j++){
  440. if(clustersPosition[j] < 0.1) continue;
  441. if(polarity*clustersCharge[j] < kClusterChargeMin) continue;
  442. bool goodHit = (clustersPosition[j]>iLo && clustersPosition[j]<iHi);
  443. double x_dut = getDUTHitPosition(j);
  444.  
  445. x_trk = x_trk0;
  446.  
  447. double dx = x_dut - x_trk;
  448. hn->Fill(dx);
  449. hnn->Fill(dx);
  450.  
  451. hca->Fill(x_trk,y_trk);
  452. if(fabs(dx)<dxWindow || iPass==1) {
  453. if(goodHit) {
  454. hx->Fill(x_trk);
  455. hxdut->Fill(x_dut);
  456. hs->Fill(clustersPosition[j]);
  457. hy->Fill(y_trk);
  458. hthx->Fill(tx);
  459. hthy->Fill(ty);
  460. hw->Fill(dx);
  461. hxy->Fill(x_trk,y_trk);
  462. ht->Fill(clustersTDC+0.1,polarity*clustersCharge[j]);
  463. hzrot->Fill(y_trk,dx);
  464. hcf->Fill(x_trk,y_trk);
  465. }
  466. }
  467. }
  468. }
  469. }
  470.  
  471. //hnn->Draw();
  472.  
  473. if(iPass==4) {
  474. findCutoutRegion(hcf);
  475. }else {
  476.  
  477. float xmin=0,xmax=0,ymin=0,ymax=0,txmin=0,txmax=0,tymin=0,tymax=0;
  478.  
  479. getRange(hthx,txmin,txmax,0.05,2);
  480. getRange(hthy,tymin,tymax,0.05,2);
  481. getRange(hx,xmin,xmax,0.05,2);
  482. getRange(hy,ymin,ymax,0.2,2);
  483. xMin = xmin;
  484. xMax = xmax;
  485. yMin = ymin;
  486. yMax = ymax;
  487. txMin = txmin;
  488. txMax = txmax;
  489. tyMin = tymin;
  490. tyMax = tymax;
  491. float xlo = 0, xhi=0;
  492. if(iPass==1) getRange(hw,xlo,xhi,0.1,1);
  493. if(iPass==2) getRange(hn,xlo,xhi,0.1,1);
  494. if(iPass==3) getRange(hnn,xlo,xhi,0.1,1);
  495. double XOFF = (xhi + xlo)/2.0;
  496. xOff = xOff - XOFF;
  497. double x_ave = (xMax + xMin)/2.0;
  498. double y_ave = (yMax + yMin)/2.0;
  499. xGloOff = xGloOff - x_ave;
  500. yGloOff = yGloOff - y_ave;
  501. // Check/correct for z rotation
  502. if(iPass==3 && correctForZRotation){
  503. cout << "Checking for z-rotation" << endl;
  504. TF1* p1 = new TF1("p1","[0]+[1]*x",yMin,yMax);
  505. p1->SetParameters(0.0,0.0001);
  506. hzrot->Fit(p1,"0R");
  507. double rz = p1->GetParameter(1);
  508. cout << "=======================================================================" << endl;
  509. cout << "==> Updating z rotation angle from " << Rz << " to " << Rz-rz << " mrad" << endl;
  510. Rz = Rz - rz;
  511. delete p1;
  512. }
  513. float lo = 0, hi = 0;
  514. getTDCBins(ht,lo,hi);
  515. tdcLo = lo;
  516. tdcHi = hi;
  517. cout << "=================================================================" << endl;
  518. cout << "====> TDC Range updated to be from " << tdcLo << " -- " << tdcHi << std::endl;
  519. cout << "=================================================================" << endl;
  520. cout << "====> Fiducial regions, xLo, xHi (Strip numbers) = " << iLo << " " << iHi << endl;
  521. cout << "====> Fiducial regions, xLo, xHi (pos in mm) = " << xMin << " " << xMax << " mm " << endl;
  522. cout << "====> Fiducial regions, yLo, yHi (pos in mm) = " << yMin << " " << yMax << " mm " << endl;
  523. cout << "====> Fiducial regions, thxLo, thxHi (pos in mrad) = " << txMin << " " << txMax << " mrad" << endl;
  524. cout << "====> Fiducial regions, thyLo, thxHi (pos in mrad) = " << tyMin << " " << tyMax << " mrad" << endl;
  525. cout << "=================================================================" << endl;
  526. cout << "====> Shifting sensor by dX = " << XOFF << " mm " << ", new xOff = " << xOff << endl;
  527. cout << "====> Global X, Y shifts: " << x_ave << " " << y_ave << endl;
  528. cout << "====> New global x,y = " << xGloOff << " " << yGloOff << endl;
  529. }
  530. //if(iPass==1) hw->Draw(); // for debugging
  531. //if(iPass==2) hy->Draw();
  532.  
  533. //return;
  534.  
  535. delete hxdut;
  536. delete hthx;
  537. delete hthy;
  538. delete hx;
  539. delete hy;
  540. delete hxy;
  541. delete hw;
  542. delete hn;
  543. delete hnn;
  544. delete hs;
  545. delete ht;
  546. delete hzrot;
  547. delete hcf;
  548. delete hca;
  549. return;
  550.  
  551. }
  552.  
  553.  
  554. Double_t AnalysisBase::getCorrChannel(double ch){
  555. //if(ch>=128 && ch<=255) return (ch + stripGap[0]);
  556. //if(ch>=256 && ch<=383) return (ch + stripGap[1]);
  557. //if(ch>=384 && ch<=512) return (ch + stripGap[2]);
  558.  
  559. if(ch<=128) return (ch - stripGap[0]);
  560. if(ch>=128 && ch<=255) return (ch - stripGap[1]);
  561. if(ch>=256 && ch<=383) return (ch - stripGap[2]);
  562. if(ch>=384 && ch<=512) return (ch - stripGap[3]);
  563.  
  564. return ch;
  565. }
  566.  
  567. void AnalysisBase::findChipBoundary(){
  568. // Look for holes in hit profile, due to boundary between chips
  569. std::cout << "=======================================================================================" << std::endl;
  570. std::cout << "findChipBoundary(): Looking for gaps due to unconnected strips in between Beetle chips " << std::endl;
  571. std::cout << "=======================================================================================" << std::endl;
  572. TH1F* hx2 = new TH1F("hx2","X position of cluster",200,-10.0,10.0);
  573.  
  574. Long64_t nbytes = 0, nb = 0;
  575. Long64_t nentries = fChain->GetEntriesFast();
  576. for (Long64_t jentry=0; jentry<max((int)nentries,20000);jentry++) {
  577. Long64_t ientry = LoadTree(jentry);
  578. if (ientry < 0) break;
  579. nb = fChain->GetEntry(jentry); nbytes += nb;
  580. for(int j=0; j<min((int)clusterNumberPerEvent,10); j++){
  581. if(clustersPosition[j] < 0.1) continue;
  582. if(polarity*clustersCharge[j] < kClusterChargeMin) continue;
  583. bool goodHit = (clustersPosition[j]>iLo && clustersPosition[j]<iHi);
  584. goodHit = goodHit && (clustersTDC >= tdcLo && clustersTDC <= tdcHi);
  585.  
  586. double x_dut = getDUTHitPosition(j);
  587. hx2->Fill(x_dut);
  588. }
  589. }
  590. //hx2->Draw(); // for debugging only
  591.  
  592. int i1 = hx2->FindBin(xMin);
  593. int i2 = hx2->FindBin(xMax);
  594. double ave = 0.0;
  595. int ixL=i1, ixH=i2;
  596. for(int i=i1; i<i2; i++){
  597. ave = ave + hx2->GetBinContent(i);
  598. }
  599. ave = ave / (i2 - i1 + 1);
  600. for(int i=i1; i<i2; i++){
  601. double v = hx2->GetBinContent(i);
  602. double xx = hx2->GetBinCenter(i);
  603. double r = v/ave;
  604. if(r < ratioForHole && xx < xLeftHole) {
  605. xLeftHole = xx;
  606. ixL = i;
  607. }
  608. if(r < ratioForHole && xx > xRightHole) {
  609. xRightHole = xx;
  610. ixH = i;
  611. }
  612. }
  613. // Count the number of "dead" channels in this region
  614. float ndead = 0;
  615. for(int i=ixL; i<=ixH; i++){
  616. double v = hx2->GetBinContent(i);
  617. double r = v/ave;
  618. if(r < ratioForHole) ndead = ndead + 1;
  619. }
  620. double frdead = ndead / (ixH - ixL + 1);
  621.  
  622. if(frdead>0.8 && ndead>4){
  623. xLeftHole = xLeftHole - 0.2;
  624. xRightHole = xRightHole + 0.2;
  625. }else{
  626. xLeftHole = 999;
  627. xRightHole = -999;
  628. }
  629. hx2->Draw();
  630. //std::cout << "Xmin, xMax = " << xMin << " " << xMax << std::endl;
  631. //std::cout << "Bin Edges: " << i1 << " " << i2 << std::endl;
  632. //std::cout << "Hole position: " << xLeftHole << " " << xRightHole << " , fraction dead = " << ndead << " " << frdead << std::endl;
  633. delete hx2;
  634.  
  635. }
  636.  
  637. void AnalysisBase::findDeadRegions(){
  638.  
  639. cout << "=================================================" << endl;
  640. cout << "findDeadRegions(): Looking for dead strip regions " << endl;
  641. cout << "=================================================" << endl;
  642.  
  643. TH1F* hf = new TH1F("hf","X position of matched cluster",200,-10.0,10.0);
  644. TH1F* hnf = new TH1F("hnf","X position of matched cluster",200,-10.0,10.0);
  645. double nomStrip = 0, detStrip = 0;
  646. Long64_t nbytes = 0, nb = 0;
  647. Int_t nentries = fChain->GetEntriesFast();
  648. for (Long64_t jentry=0; jentry<nentries;jentry++) {
  649. Long64_t ientry = LoadTree(jentry);
  650. if (ientry < 0) break;
  651. nb = fChain->GetEntry(jentry); nbytes += nb;
  652. if(n_tp3_tracks > 1) continue;
  653.  
  654. // Loop over TPIX tracks in event
  655. for(int k=0; k<n_tp3_tracks; k++){
  656. double x_trk = vec_trk_tx->at(k)*z_DUT+vec_trk_x->at(k);
  657. double y_trk = vec_trk_ty->at(k)*z_DUT+vec_trk_y->at(k);
  658.  
  659. transformTrackToDUTFrame(k, x_trk, y_trk, nomStrip, detStrip );
  660.  
  661. bool goodTrack = false;
  662. bool inFiducial = false;
  663. if(x_trk>xMin && x_trk<xMax && y_trk>yMin && y_trk<yMax) inFiducial = true;
  664. inFiducial = inFiducial && (x_trk<xLeftHole || x_trk>xRightHole);
  665. double tx = 1000*vec_trk_tx->at(k);
  666. double ty = 1000*vec_trk_ty->at(k);
  667. if(tx>txMin && tx<txMax && ty>tyMin && ty<tyMax) goodTrack = true;
  668. bool goodTime = (clustersTDC >= tdcLo && clustersTDC <= tdcHi);
  669.  
  670. bool foundHit = false;
  671. // Loop over clusters
  672. for(int j=0; j<min(clusterNumberPerEvent,10); j++){
  673. if(clustersPosition[j] < 0.1) continue;
  674. if(polarity*clustersCharge[j] < kClusterChargeMin) continue;
  675. //bool goodHit = (clustersPosition[j]>iLo || clustersPosition[j]<iHi);
  676. double x_dut = getDUTHitPosition(j);
  677.  
  678. double dx = x_dut - x_trk;
  679.  
  680. if(goodTrack && inFiducial && goodTime && fabs(dx)<dxWin) {
  681. foundHit = true;
  682. }
  683. }
  684. if(inFiducial && goodTrack && goodTime) {
  685. hf->Fill(x_trk);
  686. if(!foundHit) hnf->Fill(x_trk);
  687. }
  688. }
  689. }
  690.  
  691. TH1F *hineff = (TH1F*)hnf->Clone("hineff");
  692. hineff->Divide(hf);
  693.  
  694. // Only for debugging
  695. /*
  696. TCanvas *cc = new TCanvas("cc","Hits",800,800);
  697. cc->Divide(2,2);
  698. cc->cd(1);
  699. hnf->Draw();
  700. cc->cd(2);
  701. hf->Draw();
  702. cc->cd(3);
  703. hineff->Draw();
  704. */
  705.  
  706.  
  707. for(int k = 0; k<hineff->GetNbinsX()-1; k++){
  708. if(nDeadRegion > 20){
  709. cout << "WARNING: Exceeded 20 dead regions, exiting from loop" << endl;
  710. break;
  711. }
  712. double v1 = hineff->GetBinContent(k);
  713. if(v1 > k_DeadStrip){
  714. deadRegionLo[nDeadRegion] = hf->GetBinLowEdge(k);
  715. for(int j = k; j<hineff->GetNbinsX(); j++){
  716. double v2 = hineff->GetBinContent(j);
  717. if(v2 < k_DeadStrip){
  718. deadRegionHi[nDeadRegion] = hf->GetBinLowEdge(j);
  719. cout << "====> Found inefficiency strip region from x : "
  720. << deadRegionLo[nDeadRegion] << " <===> " << deadRegionHi[nDeadRegion] << " mm" << endl;
  721. k = j + 1;
  722. nDeadRegion++;
  723. break;
  724. }
  725. }
  726. }
  727. }
  728. delete hf;
  729. delete hnf;
  730. delete hineff;
  731. return;
  732. }
  733.  
  734. void AnalysisBase::transformTrackToDUTFrame(int k, double& x_trk, double& y_trk, double& elecStrip, double& detStrip){
  735. y_trk = y_trk + yGloOff;
  736.  
  737. double dzs = x_trk * sin(Ry);
  738. x_trk = x_trk - Rz*y_trk;
  739. y_trk = y_trk + Rz*x_trk;
  740. x_trk = x_trk + vec_trk_tx->at(k)*dzs;
  741. x_trk = x_trk/cos(Ry);
  742. double rStrip = nChan - (x_trk-xOff)/stripPitch;
  743. double corStrip = getCorrChannel(rStrip);
  744. //double dx = (corStrip - rStrip)*stripPitch;
  745. //x_trk = x_trk + dx;
  746. detStrip = rStrip; //nChan - (x_trk-xOff)/stripPitch;
  747. elecStrip = corStrip;
  748. x_trk = x_trk + xGloOff;
  749.  
  750. return;
  751. }
  752.  
  753. void AnalysisBase::PrepareDUT(){
  754.  
  755. //----------------------------
  756. // Get Beam Location (Strips)
  757. //----------------------------
  758. getBeamLoc();
  759. //return;
  760. // ==================================
  761. // Determining Optimal TDC time range
  762. // ==================================
  763. getTDC();
  764. // --------------------------------------------------------------
  765. // Find beam fiducial region (slopes & Y range) & align residuals
  766. // ---------------------------------------------------------------
  767. findBeamRegionAndAlign(1);
  768. findBeamRegionAndAlign(2);
  769. yMid = yMin + (yMax-yMin)/2.0;
  770. yHi2 = yMax - 2.0;
  771.  
  772. findBeamRegionAndAlign(3);
  773. yMid = yMin + (yMax-yMin)/2.0;
  774. yHi2 = yMax - 2.0;
  775.  
  776. //return;
  777. if(yInt2[1] > yMax){
  778. yInt1[0] = yMax-0.8; yInt1[1] = yMax-0.4;
  779. yInt2[0] = yMax-0.4; yInt2[1] = yMax;
  780. yInt3[0] = yMin; yInt3[1] = yMax-0.8;
  781. }
  782. cout << "Y Region Definitions:" << endl;
  783. cout << " ===> Region 1: " << yInt1[0] << " -- " << yInt1[1] << " mm" << endl;
  784. cout << " ===> Region 2: " << yInt2[0] << " -- " << yInt2[1] << " mm" << endl;
  785. cout << " ===> Region 3: " << yInt3[0] << " -- " << yInt3[1] << " mm" << endl;
  786. //return;
  787. // Correct for gaps between Beetle chips
  788. correctForStripGaps();
  789.  
  790. //findBeamRegionAndAlign(4);
  791. if(holeSector) findBeamRegionAndAlign(4);
  792.  
  793. findChipBoundary();
  794. cout << "====> Hole position: " << xLeftHole << " " << xRightHole << endl;
  795. //return;
  796.  
  797. if(vetoDeadStripRegions) findDeadRegions();
  798. //return;
  799.  
  800. setCrossTalkCorr();
  801.  
  802. }
  803.  
  804. Double_t AnalysisBase::getDUTHitPosition(int j){
  805. double pos = getCorrChannel(clustersPosition[j]);
  806. //double pos = clustersPosition[j];
  807. double x_dut = (nChan - pos)*stripPitch + xOff;
  808. x_dut = x_dut + xGloOff;
  809. return x_dut;
  810. }
  811.  
  812. void AnalysisBase::setCrossTalkCorr(){
  813. float biasVal = atof(m_bias);
  814. chargeCorrSlopeOdd = 0.0;
  815. chargeCorrSlopeEven = 0.0;
  816. cout << "Bias Value = " << m_bias << " " << biasVal << endl;
  817. if(m_board.Contains("A2") && m_sector=="1" && biasVal < 260){
  818. chargeCorrSlopeOdd = 0.435;
  819. chargeCorrSlopeEven = 0.400;
  820. }else if(m_board.Contains("A2") && m_sector=="1" && biasVal > 320){
  821. chargeCorrSlopeOdd = 0.148;
  822. chargeCorrSlopeEven = 0.094;
  823. }else if(m_board.Contains("A2") && m_sector=="1" && biasVal == 300){
  824. chargeCorrSlopeOdd = 0.130;
  825. chargeCorrSlopeEven = 0.088;
  826. }else if(m_board.Contains("A2") && m_sector=="2"){
  827. chargeCorrSlopeOdd = 0.120;
  828. chargeCorrSlopeEven = 0.085;
  829. }
  830. }
  831.  
  832.  
  833. void AnalysisBase::findCutoutRegion(TH2F* h){
  834.  
  835. cout << "===============================================================" << endl;
  836. cout << "findCutoutRegion(): Looking for cutout region in D type sensor " << endl;
  837. cout << "===============================================================" << endl;
  838.  
  839.  
  840. int ilx = h->GetXaxis()->FindBin(xMin) + 1;
  841. int ihx = h->GetXaxis()->FindBin(xMax) - 1;
  842. int ily = 1;//h->GetYaxis()->GetBinFindBin(yMin) - 10;
  843. int ihy = h->GetYaxis()->FindBin(yMax) + 10;
  844. if(ily < 0) ily = 0;
  845. if(ihy > h->GetYaxis()->GetNbins()) ihy = h->GetYaxis()->GetNbins();
  846. int nb = h->GetXaxis()->GetNbins();
  847.  
  848. double xlow = h->GetXaxis()->GetBinLowEdge(1);
  849. double xhi = h->GetXaxis()->GetBinLowEdge(nb)+h->GetXaxis()->GetBinWidth(1);
  850.  
  851. TH1F *hpr = new TH1F("hpr","Y vs X, edge",nb,xlow,xhi);
  852.  
  853. TH1D *hpy;
  854. int ipeak = 0;
  855. for(int i=ilx; i<ihx-1; i++){
  856. hpy = h->ProjectionY("hpy",i,i);
  857. double maxcon = 0.0;
  858. double maxbin = 0.0;
  859. for(int j=ily;j<ihy;j++){
  860. maxcon = maxcon + hpy->GetBinContent(j);
  861. if(hpy->GetBinContent(j) > maxbin) {
  862. maxbin = hpy->GetBinContent(j) ;
  863. ipeak = j;
  864. }
  865. }
  866. // Find the lower "edge"
  867. for(int j=ipeak; j>=ily; j--){
  868. double r0 = (hpy->GetBinContent(j-1)+hpy->GetBinContent(j-2))/2.0;
  869. double r1 = hpy->GetBinContent(j);
  870. double r2 = (hpy->GetBinContent(j+3)+hpy->GetBinContent(j+4))/2.0;
  871. if(r1<=2 && r0<=1.5 && r2>5*r1 && r2>8){
  872. hpr->SetBinContent(i,hpy->GetBinCenter(j)+2*hpy->GetBinWidth(j));
  873. hpr->SetBinError(i,hpy->GetBinWidth(1)/2.0);
  874. //cout << "found: " << i << " " << hpy->GetEntries() << " " << hpy->GetBinCenter(j) << " "
  875. // << hpy->GetBinContent(j) << " " << r0 << " " << r1 << " " << r2 << " " << endl;
  876. break;
  877. }
  878. }
  879. }
  880.  
  881. TF1* poly2 = new TF1("poly2","[0]+[1]*x+[2]*x*x",xMin,xMax);
  882. poly2->SetParameters(-1.5,-0.17,-0.15);
  883. hpr->Fit(poly2,"R0");
  884. hpr->SetLineColor(kRed);
  885.  
  886. holeQuadPar[0] = poly2->GetParameter(0);
  887. holeQuadPar[1] = poly2->GetParameter(1);
  888. holeQuadPar[2] = poly2->GetParameter(2);
  889.  
  890. //return;
  891.  
  892. delete poly2;
  893. delete hpr;
  894. delete hpy;
  895. return;
  896. }
  897.  
  898. bool AnalysisBase::isInCutoutRegion(double xtrk, double ytrk){
  899. // Some protections here..
  900. if(!removeTracksInHole) return false;
  901. if(holeQuadPar[0]==0 || holeQuadPar[1]==0 || holeQuadPar[2]==0) return false;
  902. // Ok, looks like we mean to really remove these tracks
  903. double yhole = holeQuadPar[0]+holeQuadPar[1]*xtrk+holeQuadPar[2]*xtrk*xtrk;
  904. if(ytrk < yhole + minDistFromHole) return true;
  905. return false;
  906. }
  907.  
  908. double AnalysisBase::DistToCutoutRegion(double xtrk, double ytrk){
  909. // Some protections here..
  910. if(holeQuadPar[0]==0 || holeQuadPar[1]==0 || holeQuadPar[2]==0) return 999.0;
  911. // Ok, looks like we mean to really remove these tracks
  912. double yhole = holeQuadPar[0]+holeQuadPar[1]*xtrk+holeQuadPar[2]*xtrk*xtrk;
  913. double a = 2.0*holeQuadPar[2]*xtrk+holeQuadPar[1];
  914. double c = yhole - a*xtrk;
  915. double b = -1.0;
  916. double dist = fabs(a*xtrk + b*ytrk + c) / sqrt(a*a+b*b);
  917. if(ytrk < yhole) dist = -1.0*dist;
  918. return dist;
  919. }
  920.  
  921. void AnalysisBase::correctForStripGaps(){
  922. cout << "================================================================================" << endl;
  923. cout << "correctForStripGaps(): Correcting x fiducial range for gaps between Beetle chips" << endl;
  924. cout << "================================================================================" << endl;
  925.  
  926. double x1 = getCorrChannel(iLo);
  927. double x2 = getCorrChannel(iHi);
  928. xMax = (nChan - x1)*stripPitch + xOff + xGloOff;
  929. xMin = (nChan - x2)*stripPitch + xOff + xGloOff;
  930. double xave = (xMax + xMin)/2.;
  931. xGloOff = xGloOff - xave;
  932. xMin = xMin - xave;
  933. xMax = xMax - xave;
  934. std::cout << "====> Updating Global Offset to use Strips, *** New xMin, xMax == " << xMin << " " << xMax << endl;
  935. std::cout << "====> xMin, xMax = " << xMin << " " << xMax << endl;
  936. std::cout << "====> xMin, xMax = " << yMin << " " << yMax << endl;
  937. return;
  938. }
  939.  
  940.