/* * File: Efficiency.cc * Author: daniel * * Created on 16. Februar 2014, 18:11 */ #include "Efficiency.h" void Efficiency::init(const TBCore* core) { for(auto dut: core->usedDUT) { int iden = dut->iden; int nCols = dut->getNcols(); int nRows = dut->getNrows(); Efficiency::ntracksraw[iden] = 0; Efficiency::ntracksgood[iden] = 0; Efficiency::ntracksregion[iden] = 0; Efficiency::ntracks[iden] = 0; Efficiency::nhits[iden] = 0; Efficiency::anyhit[iden] = 0; for (int i=0; ipixelGeometries) { double maxPitchX = geometry->getMaxPitchX(); double maxPitchY = geometry->getMaxPitchY(); double scale = maxPitchX / maxPitchY; TH2D* trackMap = new TH2D("","Track Pixel Map;Long Side [#mum];Short Side[#mum]", TMath::Nint(maxPitchX), 0, maxPitchX, TMath::Nint(maxPitchY), 0, maxPitchY); TH2D* hitMap = new TH2D("","Hit Pixel Map;Long Side [#mum];Short Side[#mum]", TMath::Nint(maxPitchX), 0, maxPitchX, TMath::Nint(maxPitchY), 0, maxPitchY); /* TH2D* trackMap = new TH2D("","Track Pixel Map;Long Side [#mum];Short Side[#mum]", TMath::Nint(maxPitchX), 0, maxPitchX, TMath::Nint(maxPitchY*scale), -maxPitchY*(scale/2), maxPitchY*(scale/2)); TH2D* hitMap = new TH2D("","Hit Pixel Map;Long Side [#mum];Short Side[#mum]", TMath::Nint(maxPitchX), 0, maxPitchX, TMath::Nint(maxPitchY*scale), -maxPitchY*(scale/2), maxPitchY*(scale/2)); */ Efficiency::h_trackPixelMap[iden].push_back(trackMap); Efficiency::h_hitPixelMap[iden].push_back(hitMap); } // get params std::string param = core->getParam(iden, Efficiency::name, "effPlotMinY"); if(param.compare("") != 0) { Efficiency::effPlotMinY[iden] = std::stod(param); } else // default value { Efficiency::effPlotMinY[iden] = 0.8; } param = core->getParam(iden, Efficiency::name, "pixMapBinSizeX"); if(param.compare("") != 0) { Efficiency::pixMapBinSizeX[iden] = std::stoi(param); } else // default value { Efficiency::pixMapBinSizeX[iden] = 8; } param = core->getParam(iden, Efficiency::name, "pixMapBinSizeY"); if(param.compare("") != 0) { Efficiency::pixMapBinSizeY[iden] = std::stoi(param); } else // default value { Efficiency::pixMapBinSizeY[iden] = 2; } } } void Efficiency::event(const TBCore* core, const TBEvent* event) { DUT* dut = event->dut; int iden = event->iden; Efficiency::ntracksraw[iden] += event->rawTracks.size(); // Check that track has passed all cuts if(event->fTracks != kGood) { // TBLOG(kINFO, "No good tracks in event "); return; } Efficiency::ntracksgood[iden] += event->tracks.size(); for(auto tbtrack: event->tracks) { // Has this track passed a masked region? if(tbtrack->fTrackMaskedRegion == kGood) { continue; } Efficiency::ntracksregion[iden]++; // Are we in the central region on the chip? if(tbtrack->fTrackEdgeRegion == kGood) { continue; } Efficiency::ntracks[iden]++; // Fill the sensor and pixel map int trackCol = tbtrack->trackCol; int trackRow = tbtrack->trackRow; Efficiency::h_trackMap[iden]->Fill(trackCol, trackRow); // Count tracks of column Efficiency::tracksvscol[iden][trackCol]++; PixelGeometry* pixelGeometry = dut->pixelArray[trackCol][trackRow].getGeometry(); int geometryNr = dut->getGeometryNumber(pixelGeometry); double foldedX; double foldedY; dut->getFoldedXY(trackCol, trackRow, tbtrack->trackX, tbtrack->trackY, &foldedX, &foldedY); Efficiency::h_trackPixelMap[iden][geometryNr]->Fill(foldedX, foldedY); if(tbtrack->fMatchedCluster == kBad) { continue; } Efficiency::nhits[iden]++; // matchedTrackMap Efficiency::h_hitMap[iden]->Fill(trackCol, trackRow); // matchedTrackPixelMap Efficiency::h_hitPixelMap[iden][geometryNr]->Fill(foldedX, foldedY); // Count hits of column Efficiency::hitsvscol[iden][trackCol]++; } // Plot the position of *any* single hits for(auto hit: event->hits) { Efficiency::h_anyHitMap[iden]->Fill(hit->col, hit->row); } } void Efficiency::finalize(const TBCore* core) { core->output->processName = Efficiency::name; char* histoTitle = new char[500]; for(auto dut: core->usedDUT) { int iden = dut->iden; int nCols = dut->getNcols(); int nRows = dut->getNrows(); // set up cuts core->output->cuts = "only events with good tracks and good clusters"; core->output->currentPreprocessCuts = ""; // ----------- std::sprintf(histoTitle, "Any Hit Map DUT %i", iden); Efficiency::h_anyHitMap[iden]->SetTitle(histoTitle); std::sprintf(histoTitle, "any_hit_map_dut_%i", iden); Efficiency::h_anyHitMap[iden]->SetName(histoTitle); core->output->drawAndSave(Efficiency::h_anyHitMap[iden], "colz", "e"); // set up cuts core->output->cuts = "only events with good tracks through central pixel and good clusters"; core->output->currentPreprocessCuts = ""; // ----------- std::sprintf(histoTitle, "Track Map DUT %i", iden); Efficiency::h_trackMap[iden]->SetTitle(histoTitle); std::sprintf(histoTitle, "track_map_dut_%i", iden); Efficiency::h_trackMap[iden]->SetName(histoTitle); core->output->drawAndSave(Efficiency::h_trackMap[iden], "colz", "e"); for(int geometryNr = 0; geometryNr < Efficiency::h_trackPixelMap[iden].size(); geometryNr++) { std::string edge = ""; if(dut->pixelGeometries[geometryNr]->isEdgePixel() == true) { edge = "(Edge Pixel)"; } std::sprintf(histoTitle, "Track Pixel Map DUT %i Geometry %i %s", iden, geometryNr, edge.c_str()); Efficiency::h_trackPixelMap[iden][geometryNr]->SetTitle(histoTitle); std::sprintf(histoTitle, "track_pixel_map_dut_%i_geometry_%i", iden, geometryNr); Efficiency::h_trackPixelMap[iden][geometryNr]->SetName(histoTitle); core->output->drawAndSave(Efficiency::h_trackPixelMap[iden][geometryNr], "colz", "e"); } // set up cuts core->output->cuts = "only good tracks through central pixel with matched cluster"; core->output->currentPreprocessCuts = ""; // ----------- std::sprintf(histoTitle, "Hit Map DUT %i", iden); Efficiency::h_hitMap[iden]->SetTitle(histoTitle); std::sprintf(histoTitle, "hit_map_dut_%i", iden); Efficiency::h_hitMap[iden]->SetName(histoTitle); core->output->drawAndSave(Efficiency::h_hitMap[iden], "colz", "e"); for(int geometryNr = 0; geometryNr < h_hitPixelMap[iden].size(); geometryNr++) { std::string edge = ""; if(dut->pixelGeometries[geometryNr]->isEdgePixel() == true) { edge = "(Edge Pixel)"; } std::sprintf(histoTitle, "Hit Pixel Map DUT %i Geometry %i %s", iden, geometryNr, edge.c_str()); Efficiency::h_hitPixelMap[iden][geometryNr]->SetTitle(histoTitle); std::sprintf(histoTitle, "hit_pixel_map_dut_%i_geometry_%i", iden, geometryNr); Efficiency::h_hitPixelMap[iden][geometryNr]->SetName(histoTitle); core->output->drawAndSave(Efficiency::h_hitPixelMap[iden][geometryNr], "colz", "e"); } TBLOG(kINFO, "DUT " << iden); TBLOG(kINFO, "Tracks raw: " << Efficiency::ntracksraw[iden]); TBLOG(kINFO, "Tracks good: " << Efficiency::ntracksgood[iden]); TBLOG(kINFO, "Tracks good region: " << Efficiency::ntracksregion[iden]); TBLOG(kINFO, "Tracks: " << Efficiency::ntracks[iden]); TBLOG(kINFO, "Tracks w/ hit: " << Efficiency::nhits[iden]); double eff = double(Efficiency::nhits[iden]) / double(Efficiency::ntracks[iden]); double eeff = TMath::Sqrt(eff * (1 - eff) / Efficiency::ntracks[iden]); TBLOG(kINFO, "Efficiency (match; no edge pixels included): " << eff << " +- " << eeff); //Calculate efficiency of columns for (int i=0; iFill(i, eff); } std::sprintf(histoTitle, "Efficiency vs Column Distribution DUT %i", iden); Efficiency::h_effvscolDis[iden]->SetTitle(histoTitle); std::sprintf(histoTitle, "eff_vs_col_dis_dut_%i", iden); Efficiency::h_effvscolDis[iden]->SetName(histoTitle); core->output->drawAndSave(Efficiency::h_effvscolDis[iden]); // Get efficiency maps TH2D* h_effMap = (TH2D*) Efficiency::h_hitMap[iden]->Clone(TString::Format("effMap")); h_effMap->SetMinimum(0.8); h_effMap->SetTitle("Efficiency Map;Column;Row"); h_effMap->Divide(Efficiency::h_hitMap[iden], Efficiency::h_trackMap[iden], 1.0, 1.0, "B"); std::sprintf(histoTitle, "Efficiency Map DUT %i", iden); h_effMap->SetTitle(histoTitle); std::sprintf(histoTitle, "efficiency_map_dut_%i", iden); h_effMap->SetName(histoTitle); core->output->drawAndSave(h_effMap, "colz", "e"); for (int geometryNr = 0; geometryNr < h_trackPixelMap[iden].size(); geometryNr++) { char* histoName = new char[600]; std::string edge = ""; if(dut->pixelGeometries[geometryNr]->isEdgePixel() == true) { edge = "(Edge Pixel)"; } Efficiency::h_trackPixelMap[iden][geometryNr]->Rebin2D(Efficiency::pixMapBinSizeX[iden], Efficiency::pixMapBinSizeY[iden]); Efficiency::h_hitPixelMap[iden][geometryNr]->Rebin2D(Efficiency::pixMapBinSizeX[iden], Efficiency::pixMapBinSizeY[iden]); std::sprintf(histoTitle, "Efficiency Pixel Map Projection DUT %i Geometry %i %s", iden, geometryNr, edge.c_str()); std::sprintf(histoName, "efficiency_pixel_map_proj_dut_%i_geometry_%i", iden, geometryNr); Efficiency::createPixelProjection(core, dut->pixelGeometries[geometryNr], iden, histoTitle, histoName, Efficiency::h_hitPixelMap[iden][geometryNr], Efficiency::h_trackPixelMap[iden][geometryNr]); std::sprintf(histoTitle, "Efficiency Pixel Map Projection Overlay DUT %i Geometry %i %s", iden, geometryNr, edge.c_str()); std::sprintf(histoName, "efficiency_pixel_map_projOver_dut_%i_geometry_%i", iden, geometryNr); Efficiency::createPixelProjectionOverlay(core, dut->pixelGeometries[geometryNr], iden, histoTitle, histoName, Efficiency::h_hitPixelMap[iden][geometryNr], Efficiency::h_trackPixelMap[iden][geometryNr]); std::sprintf(histoTitle, "Efficiency Pixel Map DUT %i Geometry %i %s", iden, geometryNr, edge.c_str()); std::sprintf(histoName, "efficiency_pixel_map_dut_%i_geometry_%i", iden, geometryNr); Efficiency::createPixelEfficiencyMap(core, dut->pixelGeometries[geometryNr], iden, histoTitle, histoName, Efficiency::h_hitPixelMap[iden][geometryNr], Efficiency::h_trackPixelMap[iden][geometryNr]); delete [] histoName; } // get inefficiency maps TH2D* h_misshitMap = (TH2D*) h_trackMap[iden]->Clone(TString::Format("misshitMap")); h_misshitMap->Add(h_trackMap[iden], h_hitMap[iden], 1.0, -1.0); TH2D* h_ineffMap = (TH2D*) h_misshitMap->Clone(TString::Format("ineffMap")); h_ineffMap->SetMinimum(0.0); h_ineffMap->SetMaximum(0.5); h_ineffMap->Divide(h_misshitMap, h_trackMap[iden], 1.0, 1.0, "B"); h_ineffMap->SetTitle("Inefficiency Map;Column;Row"); h_ineffMap->SetStats(false); std::sprintf(histoTitle, "Inefficiency Map DUT %i", iden); h_ineffMap->SetTitle(histoTitle); std::sprintf(histoTitle, "ineffMap_dut_%i", iden); h_ineffMap->SetName(histoTitle); core->output->drawAndSave(h_ineffMap, "colz", "e"); TH1D* h_effDis = new TH1D("", "eff distribution", 200, 0.0, 1.0); h_effDis->SetStats(1); TH1D* h_ineffDis = new TH1D("", "ineff distribution", 200, 0.0, 1.0); h_ineffDis->SetStats(1); TH1D* h_effDis1 = new TH1D("", "eff distribution 1", 200, 0.0, 1.0); h_effDis1->SetStats(1); TH1D* h_ineffDis1 = new TH1D("", "ineff distribution 1", 200, 0.0, 1.0); h_ineffDis1->SetStats(1); for(int ybin = 1; ybin <= nRows; ybin++) { for(int xbin = 1; xbin <= nCols; xbin++) { int bin = h_effMap->FindBin(xbin, ybin); double eff = h_effMap->GetBinContent(bin); h_effDis->Fill(eff); if(eff < 1) { h_effDis1->Fill(eff); } double ineff = h_ineffMap->GetBinContent(bin); h_ineffDis->Fill(ineff); if(ineff > 0) { h_ineffDis1->Fill(ineff); } } } std::sprintf(histoTitle, "Efficiency Distribution DUT %i", iden); h_effDis->SetTitle(histoTitle); std::sprintf(histoTitle, "effDis_dut_%i", iden); h_effDis->SetName(histoTitle); core->output->drawAndSave(h_effDis); std::sprintf(histoTitle, "Inefficiency Distribution DUT %i", iden); h_ineffDis->SetTitle(histoTitle); std::sprintf(histoTitle, "ineffDis_dut_%i", iden); h_ineffDis->SetName(histoTitle); core->output->drawAndSave(h_ineffDis); std::sprintf(histoTitle, "Efficiency Distribution 1 DUT %i", iden); h_effDis1->SetTitle(histoTitle); std::sprintf(histoTitle, "effDis1_dut_%i", iden); h_effDis1->SetName(histoTitle); core->output->drawAndSave(h_effDis1); std::sprintf(histoTitle, "Inefficiency Distribution 1 DUT %i", iden); h_ineffDis1->SetTitle(histoTitle); std::sprintf(histoTitle, "ineffDis1_dut_%i", iden); h_ineffDis1->SetName(histoTitle); core->output->drawAndSave(h_ineffDis1); delete h_effMap; delete h_misshitMap; delete h_ineffMap; delete h_effDis; delete h_ineffDis; delete h_effDis1; delete h_ineffDis1; delete Efficiency::h_effvscolDis[iden]; delete Efficiency::h_trackMap[iden]; delete Efficiency::h_hitMap[iden]; delete Efficiency::h_anyHitMap[iden]; for(auto ptr: Efficiency::h_trackPixelMap[iden]) { delete ptr; } Efficiency::h_trackPixelMap[iden].clear(); for(auto ptr: Efficiency::h_hitPixelMap[iden]) { delete ptr; } Efficiency::h_hitPixelMap[iden].clear(); delete Efficiency::h_trackPixelBiasMap[iden]; delete Efficiency::h_hitPixelBiasMap[iden]; delete Efficiency::h_trackPixelReadoutMap[iden]; delete Efficiency::h_hitPixelReadoutMap[iden]; } Efficiency::ntracksraw.clear(); Efficiency::ntracksgood.clear(); Efficiency::ntracksregion.clear(); Efficiency::ntracks.clear(); Efficiency::nhits.clear(); Efficiency::anyhit.clear(); Efficiency::tracksvscol.clear(); Efficiency::hitsvscol.clear(); Efficiency::h_effvscolDis.clear(); Efficiency::h_trackMap.clear(); Efficiency::h_hitMap.clear(); Efficiency::h_anyHitMap.clear(); Efficiency::h_trackPixelMap.clear(); Efficiency::h_hitPixelMap.clear(); Efficiency::h_trackPixelBiasMap.clear(); Efficiency::h_hitPixelBiasMap.clear(); Efficiency::h_trackPixelReadoutMap.clear(); Efficiency::h_hitPixelReadoutMap.clear(); Efficiency::effPlotMinY.clear(); Efficiency::pixMapBinSizeX.clear(); Efficiency::pixMapBinSizeY.clear(); delete[] histoTitle; } void Efficiency::createPixelProjection(const TBCore* core, const PixelGeometry* geometry, const int iden, const char* histoTitle, const char* histoName, TH2D* h_hitPixelMap, TH2D* h_trackPixelMap) { // Make projections TH1D* h_hitProj = (TH1D*) h_hitPixelMap->ProjectionX(TString::Format("%s_prjhit", h_hitPixelMap->GetName()), 1, h_hitPixelMap->GetNbinsY(), "e"); TH1D* h_trkProj = (TH1D*) h_trackPixelMap->ProjectionX(TString::Format("%s_prjtrack", h_trackPixelMap->GetName()), 1, h_trackPixelMap->GetNbinsY(), "e"); // Get the efficiency TGraphAsymmErrors* tg_effProj = new TGraphAsymmErrors(); tg_effProj->BayesDivide(h_hitProj, h_trkProj); for(int i = 0; i != tg_effProj->GetN(); i++) { tg_effProj->SetPointError(i, 0.00001, 0.000001, tg_effProj->GetErrorYlow(i), tg_effProj->GetErrorYhigh(i)); //tg_effProj->SetPointError(i,0.00001,0.000001,0.0000001,0.000001); // No errors in Y? } tg_effProj->GetXaxis()->SetRangeUser(h_hitPixelMap->GetXaxis()->GetXmin(), h_hitPixelMap->GetXaxis()->GetXmax()); // Make canvas with correct aspect ratio /*double k = geometry->getMaxPitchY() / geometry->getMaxPitchX(); TCanvas* c_proj = new TCanvas("c_proj", "", 1200, (int) floor(1200 * k)); c_proj->cd(); */ TH1D* h = (TH1D*) tg_effProj->GetHistogram(); h->GetYaxis()->SetRangeUser(Efficiency::effPlotMinY[iden], 1.01); tg_effProj->SetTitle(histoTitle); tg_effProj->SetName(histoName); std::vector hists = std::vector(); hists.push_back(new TBOutputTNamed(h, "", "", false)); hists.push_back(new TBOutputTNamed(tg_effProj, "LP", "", false)); core->output->drawAndSave(hists, histoName); delete h_hitProj; delete h_trkProj; delete h; } void Efficiency::createPixelProjectionOverlay(const TBCore* core, const PixelGeometry* geometry, const int iden, const char* histoTitle, const char* histoName, TH2D* h_hitPixelMap, TH2D* h_trackPixelMap) { // Make projections TH1D* h_hitProj = (TH1D*) h_hitPixelMap->ProjectionX(TString::Format("%s_prjhit", h_hitPixelMap->GetName()), 1, h_hitPixelMap->GetNbinsY(), "e"); TH1D* h_trkProj = (TH1D*) h_trackPixelMap->ProjectionX(TString::Format("%s_prjtrack", h_trackPixelMap->GetName()), 1, h_trackPixelMap->GetNbinsY(), "e"); h_trkProj->SetLineColor(kRed); h_trkProj->SetTitle(histoTitle); h_trkProj->SetName(histoName); std::vector hists = std::vector(); hists.push_back(new TBOutputTNamed(h_trkProj, "histo", "em", false)); hists.push_back(new TBOutputTNamed(h_hitProj, "histo sames", "em", false)); core->output->drawAndSave(hists, histoName); delete h_hitProj; delete h_trkProj; } void Efficiency::createPixelEfficiencyMap(const TBCore* core, const PixelGeometry* geometry, const int iden, const char* histoTitle, const char* histoName, TH2D* h_hitPixelMap, TH2D* h_trackPixelMap) { double k = geometry->getMaxPitchY() / geometry->getMaxPitchX(); // Get the pixel efficiency map TH2D* effPixelMap = (TH2D*) h_hitPixelMap->Clone(TString::Format("effPixelMap")); effPixelMap->Divide(h_hitPixelMap, h_trackPixelMap, 1.0, 1.0, "B"); effPixelMap->GetZaxis()->SetRangeUser(Efficiency::effPlotMinY[iden], 1.0); effPixelMap->SetTitle(histoTitle); effPixelMap->SetName(histoName); gStyle->SetPalette(1, 0); //Efficiency::SetAspectStyle(effPixelMap); /* TCanvas* canv = core->output->getTBOutputCanvas(histoName); TPad* histoPad = core->output->getTBOutputCanvasHistPad(canv); histoPad->cd(); effPixelMap->Draw("cont4z"); //effPixelMap->GetYaxis()->SetRangeUser(-30, 80); gPad->Update(); TPaletteAxis* palette = (TPaletteAxis*) effPixelMap->GetListOfFunctions()->FindObject("palette"); TPad* hist = new TPad("hist", "hist", 0, (1-k)/2, 0.87, (1-k)/2+k); if(palette != NULL) { TPad* palettePad = new TPad("palette", "palette", 0.88, 0, 1, 1); TPaletteAxis* paletteClone = (TPaletteAxis*) palette->Clone("palette"); paletteClone->SetX1NDC(0); paletteClone->SetX2NDC(0.4); paletteClone->SetY1NDC(0.03); paletteClone->SetY2NDC(0.97); paletteClone->SetLabelSize(0.2); //histoPad->Clear(); //palettePad->Draw(); //palettePad->cd(); //paletteClone->Draw(); } else { histoPad->Clear(); } histoPad->cd(); //hist->Draw(); //hist->cd(); //effPixelMap->Draw("cont4"); core->output->drawAndSaveCanvas(canv); */ core->output->drawAndSave(effPixelMap, "cont4z", "e"); delete effPixelMap; } void Efficiency::SetAspectStyle(TH1* h) { h->GetYaxis()->SetTickLength(0.01); int nDivX = (int) (h->GetXaxis()->GetXmax() / 99.0); h->GetXaxis()->SetNdivisions(nDivX); h->GetYaxis()->SetNdivisions(3); h->GetXaxis()->SetLabelOffset(0.005); h->GetYaxis()->SetLabelOffset(0.01); h->GetXaxis()->SetLabelSize(0.12); h->GetYaxis()->SetLabelSize(0.12); //h->SetTitle(";;"); return; }