/* * File: Residuals.cc * Author: daniel * * Created on 7. März 2014, 18:12 */ #include "Residuals.h" const TString Residuals::algName[] = { "eta", "qmean", "geomean", "maxtot", "DigitalHeadTail" }; void Residuals::init(const TBCore* core) { for(auto dut: core->usedDUT) { int iden = dut->iden; // get options or set standard values std::string param = core->getParam(iden, Residuals::name, "use_minimum_sum_of_tot_per_cluster"); if(param.compare("false") == 0) { Residuals::use_minimum_sum_of_tot_per_cluster[iden] = false; } else // default value { Residuals::use_minimum_sum_of_tot_per_cluster[iden] = true; } param = core->getParam(iden, Residuals::name, "minimum_sum_of_tot_per_cluster"); if(param.compare("") != 0) { Residuals::minimum_sum_of_tot_per_cluster[iden] = std::stoi(param); } else //default value { Residuals::minimum_sum_of_tot_per_cluster[iden] = 3; } param = core->getParam(iden, Residuals::name, "use_minimum_entries_to_plot_histo"); if(param.compare("true") == 0) { Residuals::use_minimum_entries_to_plot_histo[iden] = true; } else // default value { Residuals::use_minimum_entries_to_plot_histo[dut->iden] = false; } param = core->getParam(iden, Residuals::name, "minimum_entries_to_plot_histo"); if(param.compare("") != 0) { Residuals::minimum_entries_to_plot_histo[dut->iden] = std::stoi(param); } else // default value { Residuals::minimum_entries_to_plot_histo[dut->iden] = 20; } int nCols = dut->getNcols(); int nRows = dut->getNrows(); double maxPixelPitchX = dut->getMaxPixelPitchX(); double maxPixelPitchY = dut->getMaxPixelPitchY(); // initialize standard histograms (all other histograms are created dynamically as need arises) Residuals::h_angleX[iden] = new TH1D("", ";Beam dx/dz;Entries", 1000, -0.001, 0.001); Residuals::h_angleY[iden] = new TH1D("", ";Beam dy/dz;Entries", 1000, -0.001, 0.001); Residuals::h_clusizeX[iden] = new TH1D("", ";Cluster span [col];Entries", nCols, -.5, nCols-.5); Residuals::h_clusizeY[iden] = new TH1D("", ";Cluster span [row];Entries", nRows, -.5, nRows-.5); Residuals::h_clusize[iden] = new TH1D("", ";Cluster size;Entries", nRows, -.5, nRows-.5); Residuals::h_clusizeXY[iden] = new TH2D("", ";Cluster span [col];Cluster span [row]", nCols, -.5, nCols-.5, nRows, -.5, nRows-.5); // fuer einzelne Pixel evtl ??? //Residuals::h_ecorrsX[iden] = new TH1D("", ";Eta Correction X [#mum];Entries", 100, 0, maxPixelPitchX); //Residuals::h_ecorrsY[iden] = new TH1D("", ";Eta Correction Y [#mum];Entries", 100, 0, maxPixelPitchY); } } void Residuals::event(const TBCore* core, const TBEvent* event) { DUT* dut = event->dut; int iden = event->iden; if(event->fTracks != kGood) { return; } for(auto track: event->tracks) { Residuals::h_angleX[iden]->Fill(track->dxdz); Residuals::h_angleY[iden]->Fill(track->dydz); if(track->fTrackMaskedRegion == kGood) { continue; } if(track->fMatchedCluster != kGood) { continue; } TBCluster* cluster = track->matchedCluster; // Check that cluster has charge if(Residuals::use_minimum_sum_of_tot_per_cluster[iden] == true) { if(TBCluster::getSumToT(cluster) < Residuals::minimum_sum_of_tot_per_cluster[iden]) { continue; } // make FE-dependent ToT cut? } // Plot residuals based on cluster size int nCols = TBCluster::spannedCols(cluster, event); int nRows = TBCluster::spannedRows(cluster, event); int clusterSize = cluster->hits.size(); Residuals::h_clusizeX[iden]->Fill(nCols); Residuals::h_clusizeY[iden]->Fill(nRows); Residuals::h_clusize[iden]->Fill(cluster->hits.size()); Residuals::h_clusizeXY[iden]->Fill(nCols, nRows); // TODO think about how to make it better double maxPixelPitchX = dut->getMaxPixelPitchX(); double maxPixelPitchY = dut->getMaxPixelPitchY(); char* histoTitle = new char[500]; // check whether there was already created a histogram for the cluster size nCols in X --> this cluster size already occurred and we can just fill the histogram if(Residuals::hhh_resX[iden][nCols].empty()) { for(int method = 0; method < algNum; method++) { std::sprintf(histoTitle, "%s_cluSize%i_resX_dut_%i;X Residual [#mum];Entries", Residuals::algName[method].Data(), nCols, iden); Residuals::hhh_resX[iden][nCols].push_back(new TH1D("", histoTitle, 100, -maxPixelPitchX * 2 * nCols, maxPixelPitchX * 2 * nCols)); // TODO Residuals::hhh_resX[iden][nCols][method]->SetDirectory(0); } } // check whether there was already created a histogram for the cluster size nRows in Y --> this cluster size already occurred and we can just fill the histogram if(Residuals::hhh_resY[iden][nRows].empty()) { for(int method = 0; method < algNum; method++) { std::sprintf(histoTitle, "cluSize%i_resY_dut_%i;Y Residual [#mum];Entries", nRows, iden); Residuals::hhh_resY[iden][nRows].push_back(new TH1D("", histoTitle, 100, -(maxPixelPitchY * 2 * nRows), maxPixelPitchY * 2 * nRows)); // TODO Residuals::hhh_resY[iden][nRows][method]->SetDirectory(0); } } // check whether there was already created a histogram for the total cluster size cluster.size() --> this cluster size already occurred and we can just fill the histogram if(Residuals::hhh_res[iden][clusterSize].empty()) { for(int method = 0; method < algNum; method++) { std::sprintf(histoTitle, "cluSize%i_res_dut_%i;Residual [#mum];Entries", (int) cluster->hits.size(), iden); Residuals::hhh_res[iden][clusterSize].push_back(new TH1D("", histoTitle, 100, -(max(2 * maxPixelPitchY * nRows, maxPixelPitchX * nCols)), max(2 * maxPixelPitchY * nRows, maxPixelPitchX * nCols))); // TODO Residuals::hhh_res[iden][clusterSize][method]->SetDirectory(0); } } Residuals::getResiduals(event, track, &Residuals::hhh_resX[iden][nCols], &Residuals::hhh_resY[iden][nRows], &Residuals::hhh_res[iden][clusterSize]); /** * ETA correction */ delete [] histoTitle; } } void Residuals::finalize(const TBCore* core) { core->output->processName = Residuals::name; char* histoTitle = new char[500]; char* histoName = new char[500]; for(auto dut: core->usedDUT) { int iden = dut->iden; double maxPixelPitchX = dut->getMaxPixelPitchX(); double maxPixelPitchY = dut->getMaxPixelPitchY(); // set up cuts if(Residuals::use_minimum_sum_of_tot_per_cluster[iden] == true) { core->output->cuts = "only matched tracks with ToT sum >= "+std::to_string(Residuals::minimum_sum_of_tot_per_cluster[iden]); } else { core->output->cuts = "only matched tracks with matched cluster size <= 4"; } core->output->currentPreprocessCuts = ""; // ----------- std::sprintf(histoTitle, "Angle X DUT %i", iden); Residuals::h_angleX[iden]->SetTitle(histoTitle); std::sprintf(histoTitle, "angleX_dut_%i", iden); Residuals::h_angleX[iden]->SetName(histoTitle); core->output->drawAndSave(Residuals::h_angleX[iden], "", "e"); std::sprintf(histoTitle, "Angle Y DUT %i", iden); Residuals::h_angleY[iden]->SetTitle(histoTitle); std::sprintf(histoTitle, "angleY_dut_%i", iden); Residuals::h_angleY[iden]->SetName(histoTitle); core->output->drawAndSave(Residuals::h_angleY[iden], "", "e"); std::sprintf(histoTitle, "Cluster Size X DUT %i", iden); Residuals::h_clusizeX[iden]->SetTitle(histoTitle); std::sprintf(histoTitle, "clusizeX_dut_%i", iden); Residuals::h_clusizeX[iden]->SetName(histoTitle); core->output->drawAndSave(Residuals::h_clusizeX[iden], "", "em"); std::sprintf(histoTitle, "Cluster Size Y DUT %i", iden); Residuals::h_clusizeY[iden]->SetTitle(histoTitle); std::sprintf(histoTitle, "clusizeY_dut_%i", iden); Residuals::h_clusizeY[iden]->SetName(histoTitle); core->output->drawAndSave(Residuals::h_clusizeY[iden], "", "em"); std::sprintf(histoTitle, "Cluster Size DUT %i", iden); Residuals::h_clusize[iden]->SetTitle(histoTitle); std::sprintf(histoTitle, "clusize_dut_%i", iden); Residuals::h_clusize[iden]->SetName(histoTitle); core->output->drawAndSave(Residuals::h_clusize[iden], "", "em"); std::sprintf(histoTitle, "Cluster Size XY DUT %i", iden); Residuals::h_clusizeXY[iden]->SetTitle(histoTitle); std::sprintf(histoTitle, "clusizeXY_dut_%i", iden); Residuals::h_clusizeXY[iden]->SetName(histoTitle); core->output->drawAndSave(Residuals::h_clusizeXY[iden], "", "em"); /*std::sprintf(histoTitle, "Eta Correction X DUT %i", iden); Residuals::h_ecorrsX[iden]->SetTitle(histoTitle); std::sprintf(histoTitle, "etacorrsX_dut_%i", iden); core->tboutput->drawAndSave(Residuals::name, histoTitle, Residuals::h_ecorrsX[iden]); std::sprintf(histoTitle, "Eta Correction Y DUT %i", iden); Residuals::h_ecorrsY[iden]->SetTitle(histoTitle); std::sprintf(histoTitle, "etacorrsY_dut_%i", iden); core->tboutput->drawAndSave(Residuals::name, histoTitle, Residuals::h_ecorrsY[iden]); */ for(auto resX: Residuals::hhh_resX[iden]) { int size = resX.first; for(int alg = 0; alg < Residuals::algNum; alg++) { if((Residuals::use_minimum_entries_to_plot_histo[iden] == false) or (Residuals::hhh_resX[iden][size][alg]->GetEntries() >= Residuals::minimum_entries_to_plot_histo[iden])) { std::sprintf(histoTitle, "Cluster Size X %i %s Algorithm DUT %i", size, Residuals::algName[alg].Data(), iden); Residuals::hhh_resX[iden][size][alg]->SetTitle(histoTitle); std::sprintf(histoTitle, "clusizeX_%i_with_%s_algorithm_dut_%i", size, Residuals::algName[alg].Data(), iden); Residuals::hhh_resX[iden][size][alg]->SetName(histoTitle); core->output->drawAndSave(Residuals::hhh_resX[iden][size][alg], "", "emr"); std::sprintf(histoTitle, "Fit Residual Cluster Size X %i %s Algorithm DUT %i", size, Residuals::algName[alg].Data(), iden); std::sprintf(histoName, "fitResidual_clusizeX_%i_with_%s_algorithm_dut_%i", size, Residuals::algName[alg].Data(), iden); Residuals::fitResiduals(core, histoTitle, histoName, Residuals::hhh_resX[iden][size][alg], maxPixelPitchX/2.0, -maxPixelPitchX*size, maxPixelPitchX*size); std::sprintf(histoTitle, "Fit Residual multiG Cluster Size X %i %s Algorithm DUT %i", size, Residuals::algName[alg].Data(), iden); std::sprintf(histoName, "fitResidual_multiG_clusizeX_%i_with_%s_algorithm_dut_%i", size, Residuals::algName[alg].Data(), iden); Residuals::fitAllClusResiduals(core, histoTitle, histoName, Residuals::hhh_resX[iden][size][alg], maxPixelPitchX/2.0, -maxPixelPitchX*size, maxPixelPitchX*size); } } std::sprintf(histoTitle, "Comparison Between Algorithm Cluster Size X %i DUT %i", size, iden); std::sprintf(histoName, "comparison_between_algorithm_clusizeX_%i_dut_%i", size, iden); Residuals::compareResidualsForOneClustersize(core, iden, histoTitle, histoName, &Residuals::hhh_resX[iden][size]); for(int alg = 0; alg < Residuals::algNum; alg++) { delete Residuals::hhh_resX[iden][size][alg]; } } for(auto resY: Residuals::hhh_resY[iden]) { int size = resY.first; for(int alg = 0; alg < Residuals::algNum; alg++) { if((Residuals::use_minimum_entries_to_plot_histo[iden] == false) or (Residuals::hhh_resY[iden][size][alg]->GetEntries() >= Residuals::minimum_entries_to_plot_histo[iden])) { std::sprintf(histoTitle, "Cluster Size Y %i %s Algorithm DUT %i", size, Residuals::algName[alg].Data(), iden); Residuals::hhh_resY[iden][size][alg]->SetTitle(histoTitle); std::sprintf(histoTitle, "clusizeY_%i_with_%s_algorithm_dut_%i", size, Residuals::algName[alg].Data(), iden); Residuals::hhh_resY[iden][size][alg]->SetName(histoTitle); core->output->drawAndSave(Residuals::hhh_resY[iden][size][alg], "", "emr"); std::sprintf(histoTitle, "Fit Residual Cluster Size Y %i %s Algorithm DUT %i", size, Residuals::algName[alg].Data(), iden); std::sprintf(histoName, "fitResidual_clusizeY_%i_with_%s_algorithm_dut_%i", size, Residuals::algName[alg].Data(), iden); Residuals::fitResiduals(core, histoTitle, histoName, Residuals::hhh_resY[iden][size][alg], maxPixelPitchY/2.0, -maxPixelPitchY*size, maxPixelPitchY*size); std::sprintf(histoTitle, "Fit Residual multiG Cluster Size Y %i %s Algorithm DUT %i", size, Residuals::algName[alg].Data(), iden); std::sprintf(histoName, "fitResidual_multiG_clusizeY_%i_with_%s_algorithm_dut_%i", size, Residuals::algName[alg].Data(), iden); Residuals::fitAllClusResiduals(core, histoTitle, histoName, Residuals::hhh_resY[iden][size][alg], maxPixelPitchY/2.0, -maxPixelPitchY*size, maxPixelPitchY*size); } } std::sprintf(histoTitle, "Comparison Between Algorithm Cluster Size Y %i DUT %i", size, iden); std::sprintf(histoName, "comparison_between_algorithm_clusizeY_%i_dut_%i", size, iden); Residuals::compareResidualsForOneClustersize(core, iden, histoTitle, histoName, &Residuals::hhh_resY[iden][size]); for(int alg = 0; alg < Residuals::algNum; alg++) { delete Residuals::hhh_resY[iden][size][alg]; } } for(auto res: Residuals::hhh_res[iden]) { int size = res.first; for(int alg = 0; alg < Residuals::algNum; alg++) { if((Residuals::use_minimum_entries_to_plot_histo[iden] == false) or (Residuals::hhh_res[iden][size][alg]->GetEntries() >= Residuals::minimum_entries_to_plot_histo[iden])) { std::sprintf(histoTitle, "Cluster Size %i %s Algorithm DUT %i", size, Residuals::algName[alg].Data(), iden); Residuals::hhh_res[iden][size][alg]->SetTitle(histoTitle); std::sprintf(histoTitle, "clusize_%i_with_%s_algorithm_dut_%i", size, Residuals::algName[alg].Data(), iden); Residuals::hhh_res[iden][size][alg]->SetName(histoTitle); core->output->drawAndSave(Residuals::hhh_res[iden][size][alg], "", "emr"); std::sprintf(histoTitle, "Fit Residual Cluster Size %i %s Algorithm DUT %i", size, Residuals::algName[alg].Data(), iden); std::sprintf(histoName, "fitResidual_clusize_%i_with_%s_algorithm_dut_%i", size, Residuals::algName[alg].Data(), iden); Residuals::fitResiduals(core, histoTitle, histoName, Residuals::hhh_res[iden][size][alg], std::max(maxPixelPitchX / 2.0, maxPixelPitchY / 2.0), -(std::max(maxPixelPitchX * size, 2 * maxPixelPitchY * size)), std::max(maxPixelPitchX * size, 2 * maxPixelPitchY * size)); std::sprintf(histoTitle, "Fit Residual multiG Cluster Size %i %s Algorithm DUT %i", size, Residuals::algName[alg].Data(), iden); std::sprintf(histoName, "fitResidual_multiG_clusize_%i_with_%s_algorithm_dut_%i", size, Residuals::algName[alg].Data(), iden); Residuals::fitAllClusResiduals(core, histoTitle, histoName, Residuals::hhh_res[iden][size][alg], std::max(maxPixelPitchX / 2.0, maxPixelPitchY / 2.0), -(std::max(maxPixelPitchX * size, 2 * maxPixelPitchY * size)), std::max(maxPixelPitchX * size, 2 * maxPixelPitchY * size)); } } std::sprintf(histoTitle, "Comparison Between Algorithm Cluster Size %i DUT %i", size, iden); std::sprintf(histoName, "comparison_between_algorithm_clusize_%i_dut_%i", size, iden); Residuals::compareResidualsForOneClustersize(core, iden, histoTitle, histoName, &Residuals::hhh_res[iden][size]); for(int alg = 0; alg < Residuals::algNum; alg++) { delete Residuals::hhh_res[iden][size][alg]; } } delete Residuals::h_angleX[iden]; delete Residuals::h_angleY[iden]; delete Residuals::h_clusizeX[iden]; delete Residuals::h_clusizeY[iden]; delete Residuals::h_clusize[iden]; delete Residuals::h_clusizeXY[iden]; //delete Residuals::h_ecorrsX[iden]; //delete Residuals::h_ecorrsY[iden]; } Residuals::h_angleX.clear(); Residuals::h_angleY.clear(); Residuals::h_clusizeX.clear(); Residuals::h_clusizeY.clear(); Residuals::h_clusize.clear(); Residuals::h_clusizeXY.clear(); //Residuals::h_ecorrsX.clear(); //Residuals::h_ecorrsY.clear(); delete [] histoTitle; delete [] histoName; } void Residuals::getResiduals(const TBEvent* event, const TBTrack* track, std::vector* histosX, std::vector* histosY, std::vector* histosAll) { TBCluster* cluster = track->matchedCluster; TBHit* maxCell = TBCluster::getMaxToTCell(cluster); double distance; //DHT (*histosX)[4]->Fill(track->trackX - TBCluster::getClusterCenterXByDigitalHeadTail1D(cluster, event)); (*histosY)[4]->Fill(track->trackY- TBCluster::getClusterCenterYByDigitalHeadTail1D(cluster, event)); distance = TMath::Sqrt(TMath::Power((track->trackX - TBCluster::getClusterCenterXByDigitalHeadTail1D(cluster, event)), 2) + TMath::Power(track->trackY - TBCluster::getClusterCenterYByDigitalHeadTail1D(cluster, event), 2)); (*histosAll)[4]->Fill(distance); //maxcell int col = maxCell->col; int row = maxCell->row; (*histosX)[3]->Fill(track->trackX - event->dut->pixelArray[col][row].getGeometricCenterX()); // can be negative (*histosY)[3]->Fill(track->trackY - event->dut->pixelArray[col][row].getGeometricCenterY()); // can be negative distance = TMath::Sqrt(TMath::Power((track->trackY - event->dut->pixelArray[col][row].getGeometricCenterY()), 2) + TMath::Power(track->trackX - event->dut->pixelArray[col][row].getGeometricCenterX(), 2)); (*histosAll)[3]->Fill(distance); //geomean (*histosX)[2]->Fill(track->trackX - TBCluster::getUnWeightedX(cluster, event)); (*histosY)[2]->Fill(track->trackY - TBCluster::getUnWeightedY(cluster, event)); distance = TMath::Sqrt(TMath::Power((track->trackX - TBCluster::getUnWeightedX(cluster, event)), 2) + TMath::Power(track->trackY - TBCluster::getUnWeightedY(cluster, event), 2)); (*histosAll)[2]->Fill(distance); //qmean (*histosX)[1]->Fill(track->trackX - TBCluster::getChargeWeightedX(cluster, event)); (*histosY)[1]->Fill(track->trackY - TBCluster::getChargeWeightedY(cluster, event)); distance = TMath::Sqrt(TMath::Power((track->trackX - TBCluster::getChargeWeightedX(cluster, event)), 2) + TMath::Power(track->trackY - TBCluster::getChargeWeightedY(cluster, event), 2)); (*histosAll)[1]->Fill(distance); /** * ETA correction */ } void Residuals::fitAllClusResiduals(const TBCore* core, const char* histoTitle, const char* histoName, TH1* h_residuals, double halfPitch, double xmin, double xmax) { TH1* h_residualsClone = (TH1*) h_residuals->Clone(); h_residualsClone->SetTitle(histoTitle); h_residualsClone->SetName(histoName); TF1* multig = new TF1("multiG", multiG, xmin, xmax, 6); double coreMeanGuess = 0.0; double outMeanGuess = 0.0; double coreWidthGuess = halfPitch * 2.0 / sqrt(12.); double outWidthGuess = 6.0 * (halfPitch * 2.0); double coreFracGuess = 0.7; double coreFracMin = 0.0; double coreFracMax = 1.0; double normGuess = h_residualsClone->GetEntries(); int coreMeanIndex = 0; int outMeanIndex = 1; int coreWidthIndex = 2; int outWidthIndex = 3; int coreFracIndex = 4; int normIndex = 5; multig->SetParameters(coreMeanGuess, outMeanGuess, coreWidthGuess, outWidthGuess, coreFracGuess, normGuess); multig->FixParameter(outMeanIndex, outMeanGuess); multig->SetParLimits(coreFracIndex, coreFracMin, coreFracMax); multig->SetLineColor(kRed); h_residualsClone->Fit(multig, "Q +"); TCanvas* canvas = core->output->getTBOutputCanvas(histoName); h_residualsClone->SetStats(kTRUE); gStyle->SetOptStat("emR"); TPad* hist = core->output->getTBOutputCanvasHistPad(canvas); TPad* statPad = core->output->getTBOutputCanvasStatPad(canvas); hist->cd(); h_residualsClone->Draw(); gPad->Update(); // stat box TPaveStats* st; st = (TPaveStats*) h_residualsClone->GetListOfFunctions()->FindObject("stats")->Clone("stats"); st->SetX1NDC(0); st->SetY1NDC(0.01); st->SetX2NDC(0.98); st->SetY2NDC(0.99); statPad->cd(); st->Draw(); ((TH1*) canvas->FindObject(h_residualsClone->GetName()))->SetStats(kFALSE); hist->cd(); // param box multiG TF1* fitMultiG = h_residualsClone->GetFunction("multiG"); if(fitMultiG != NULL) { TPaveStats* paramBoxFunc = new TPaveStats(0.6, 0.13, 0.89, 0.4, "NDC"); paramBoxFunc->AddText("Box Fit multiG"); char* nameParam[4] = { (char*) "core sigma", (char*) "outlier sigma", (char*) "core fraction", (char*) "outlier fraction" }; char* line = new char[300]; for(int i = 0; i < 4; i++) { // TODO not sure about last param (outlier fraction) in tbmon 1 it was calculate std::sprintf(line, "%s = %10.2e #pm %10.2e", nameParam[i], fitMultiG->GetParameter(i+2), fitMultiG->GetParError(i+2)); paramBoxFunc->AddText(line); } paramBoxFunc->SetFillColor(0); paramBoxFunc->SetTextColor(kRed); paramBoxFunc->SetBorderSize(1); paramBoxFunc->Draw(); delete [] line; } core->output->setStandardTBOutputText(canvas, ""); core->output->drawAndSaveCanvas(canvas); core->output->tOutFile->cd(Residuals::name.c_str()); h_residualsClone->Write(); delete multig; delete fitMultiG; delete h_residualsClone; return; } void Residuals::compareResidualsForOneClustersize(const TBCore* core, int iden, const char* histoTitle, const char* histoName, std::vector* hh) { double maxval = (*hh)[MAXTOT]->GetMaximum(); if((*hh)[GEOMEAN]->GetMaximum() > maxval) { maxval = (*hh)[GEOMEAN]->GetMaximum(); } if((*hh)[QMEAN]->GetMaximum() > maxval) { maxval = (*hh)[QMEAN]->GetMaximum(); } if((*hh)[ETA]->GetMaximum() > maxval) { maxval = (*hh)[ETA]->GetMaximum(); } if((*hh)[DHT]->GetMaximum() > maxval) { maxval = (*hh)[DHT]->GetMaximum(); } double lw = 2.0; TCanvas* canvas = core->output->getTBOutputCanvas(histoName); TPad* hist = core->output->getTBOutputCanvasHistPad(canvas); hist->cd(); (*hh)[MAXTOT]->SetMaximum(maxval * 1.1); (*hh)[MAXTOT]->SetStats(false); (*hh)[MAXTOT]->SetLineWidth(lw); (*hh)[MAXTOT]->SetLineColor(kOrange); (*hh)[MAXTOT]->SetFillColor(kOrange); (*hh)[MAXTOT]->SetFillStyle(1001); (*hh)[MAXTOT]->SetTitle(histoTitle); if((Residuals::use_minimum_entries_to_plot_histo[iden] == false) || ((*hh)[MAXTOT]->GetEntries() >= Residuals::minimum_entries_to_plot_histo[iden])) { (*hh)[MAXTOT]->Draw(); } (*hh)[ETA]->SetLineWidth(lw); (*hh)[ETA]->SetStats(false); (*hh)[ETA]->SetLineColor(kBlack); (*hh)[ETA]->SetFillColor(kBlack); (*hh)[ETA]->SetFillStyle(3004); (*hh)[ETA]->SetTitle(histoTitle); if((Residuals::use_minimum_entries_to_plot_histo[iden] == false) || ((*hh)[ETA]->GetEntries() >= Residuals::minimum_entries_to_plot_histo[iden])) { (*hh)[ETA]->Draw("same"); } (*hh)[GEOMEAN]->SetLineWidth(lw); (*hh)[GEOMEAN]->SetStats(true); (*hh)[GEOMEAN]->SetLineColor(kBlue); (*hh)[GEOMEAN]->SetTitle(histoTitle); if((Residuals::use_minimum_entries_to_plot_histo[iden] == false) || ((*hh)[GEOMEAN]->GetEntries() >= Residuals::minimum_entries_to_plot_histo[iden])) { (*hh)[GEOMEAN]->Draw("same"); } (*hh)[QMEAN]->SetLineWidth(lw); (*hh)[QMEAN]->SetStats(false); (*hh)[QMEAN]->SetLineColor(kRed); (*hh)[QMEAN]->SetTitle(histoTitle); if((Residuals::use_minimum_entries_to_plot_histo[iden] == false) || ((*hh)[QMEAN]->GetEntries() >= Residuals::minimum_entries_to_plot_histo[iden])) { (*hh)[QMEAN]->Draw("same"); } (*hh)[DHT]->SetLineWidth(lw); (*hh)[DHT]->SetStats(false); (*hh)[DHT]->SetLineColor(kGreen); (*hh)[DHT]->SetTitle(histoTitle); if((Residuals::use_minimum_entries_to_plot_histo[iden] == false) || ((*hh)[DHT]->GetEntries() >= Residuals::minimum_entries_to_plot_histo[iden])) { (*hh)[DHT]->Draw("same"); } TLegend* leg = new TLegend(0, 0.01, 0.98, 0.99); leg->SetFillColor(kWhite); leg->SetBorderSize(1); if((Residuals::use_minimum_entries_to_plot_histo[iden] == false) || ((*hh)[MAXTOT]->GetEntries() >= Residuals::minimum_entries_to_plot_histo[iden])) { leg->AddEntry((*hh)[MAXTOT], "Max ToT", "F"); } if((Residuals::use_minimum_entries_to_plot_histo[iden] == false) || ((*hh)[ETA]->GetEntries() >= Residuals::minimum_entries_to_plot_histo[iden])) { leg->AddEntry((*hh)[ETA], "#eta corr", "F"); } if((Residuals::use_minimum_entries_to_plot_histo[iden] == false) || ((*hh)[GEOMEAN]->GetEntries() >= Residuals::minimum_entries_to_plot_histo[iden])) { leg->AddEntry((*hh)[GEOMEAN], "Digital", "L"); } if((Residuals::use_minimum_entries_to_plot_histo[iden] == false) || ((*hh)[QMEAN]->GetEntries() >= Residuals::minimum_entries_to_plot_histo[iden])) { leg->AddEntry((*hh)[QMEAN], "Analog", "L"); } if((Residuals::use_minimum_entries_to_plot_histo[iden] == false) || ((*hh)[DHT]->GetEntries() >= Residuals::minimum_entries_to_plot_histo[iden])) { leg->AddEntry((*hh)[DHT], "DHT", "L"); } TPad* stat = core->output->getTBOutputCanvasStatPad(canvas); stat->cd(); leg->Draw(); core->output->setStandardTBOutputText(canvas, ""); if(!((Residuals::use_minimum_entries_to_plot_histo[iden] == true) && ((*hh)[ETA]->GetEntries() < Residuals::minimum_entries_to_plot_histo[iden]) && ((*hh)[GEOMEAN]->GetEntries() < Residuals::minimum_entries_to_plot_histo[iden]) && ((*hh)[QMEAN]->GetEntries() < Residuals::minimum_entries_to_plot_histo[iden]) && ((*hh)[MAXTOT]->GetEntries() < Residuals::minimum_entries_to_plot_histo[iden]))) { core->output->drawAndSaveCanvas(canvas); } delete leg; return; } void Residuals::fitResiduals(const TBCore* core, const char* histoTitle, const char* histoName, TH1* h_residuals, double halfPitch, double xmin, double xmax) { TH1* h_residualsClone = (TH1*) h_residuals->Clone(); h_residualsClone->SetTitle(histoTitle); h_residualsClone->SetName(histoName); TF1* func = new TF1("func", fitFunc, xmin, xmax, 5); func->SetParameters(8.0, h_residualsClone->Integral("width"), 0.0, halfPitch, 0.0); func->FixParameter(3, halfPitch); func->SetLineColor(2); h_residualsClone->Fit(func, "+ Q"); TF1* funcBox = new TF1("funcBox", fitFuncBox, xmin, xmax, 3); funcBox->SetParameters(8.0, h_residualsClone->Integral("width"), halfPitch); funcBox->FixParameter(3, halfPitch); funcBox->SetLineColor(kBlue); h_residualsClone->Fit(funcBox, "+ Q"); TF1* funcGaus = new TF1("funcGaus", "gaus", xmin, xmax); funcGaus->SetParameters(8.0, h_residualsClone->Integral("width"), halfPitch); funcGaus->FixParameter(3, halfPitch); funcGaus->SetLineColor(kGreen); h_residualsClone->Fit(funcGaus, "+ Q"); TCanvas* canvas = core->output->getTBOutputCanvas(histoName); h_residualsClone->SetStats(kTRUE); gStyle->SetOptStat("emR"); TPad* hist = core->output->getTBOutputCanvasHistPad(canvas); TPad* statPad = core->output->getTBOutputCanvasStatPad(canvas); hist->cd(); h_residualsClone->Draw(); gPad->Update(); // stat box TPaveStats* st; st = (TPaveStats*) h_residualsClone->GetListOfFunctions()->FindObject("stats")->Clone("stats"); st->SetX1NDC(0); st->SetY1NDC(0.01); st->SetX2NDC(0.98); st->SetY2NDC(0.99); statPad->cd(); st->Draw(); ((TH1*) canvas->FindObject(h_residualsClone->GetName()))->SetStats(kFALSE); hist->cd(); // param box func TF1* fitFunc = h_residualsClone->GetFunction("func"); if(fitFunc != NULL) { TPaveStats* paramBoxFunc = new TPaveStats(0.6, 0.13, 0.89, 0.4, "NDC"); paramBoxFunc->AddText("Box fit (w / slope)"); char* nameParam[5] = { (char*) "sigma", (char*) "area", (char*) "x0", (char*) "width/2", (char*) "angCoeff" }; char* line = new char[300]; for(int i = 0; i < 5; i++) { std::sprintf(line, "%s = %10.2e #pm %10.2e", nameParam[i], fitFunc->GetParameter(i), fitFunc->GetParError(i)); paramBoxFunc->AddText(line); } paramBoxFunc->SetFillColor(0); paramBoxFunc->SetTextColor(kRed); paramBoxFunc->SetBorderSize(1); paramBoxFunc->Draw(); delete [] line; } // param box funcBox TF1* fitFuncBox = h_residualsClone->GetFunction("funcBox"); if(fitFuncBox != NULL) { TPaveStats *paramBoxFunc = new TPaveStats(0.6, 0.42, 0.89, 0.62, "NDC"); paramBoxFunc->AddText("Box fit (w/o slope)"); char* nameParam[3] = { (char*) "sigma", (char*) "area", (char*) "width/2" }; char* line = new char[300]; for(int i = 0; i < 3; i++) { std::sprintf(line, "%s = %10.2e #pm %10.2e", nameParam[i], fitFuncBox->GetParameter(i), fitFuncBox->GetParError(i)); paramBoxFunc->AddText(line); } paramBoxFunc->SetFillColor(0); paramBoxFunc->SetTextColor(kBlue); paramBoxFunc->SetBorderSize(1); paramBoxFunc->Draw(); delete [] line; } // param box gaus TF1* fitFuncGaus = h_residualsClone->GetFunction("funcGaus"); if(fitFuncGaus != NULL) { TPaveStats *paramBoxFunc = new TPaveStats(0.6, 0.64, 0.89, 0.84, "NDC"); paramBoxFunc->AddText("Gaus Fit"); char* nameParam[3] = { (char*) "area", (char*) "mean", (char*) "sigma" }; char* line = new char[300]; for(int i = 0; i < 3; i++) { std::sprintf(line, "%s = %10.2e #pm %10.2e", nameParam[i], fitFuncGaus->GetParameter(i), fitFuncGaus->GetParError(i)); paramBoxFunc->AddText(line); } paramBoxFunc->SetFillColor(0); paramBoxFunc->SetTextColor(kGreen); paramBoxFunc->SetBorderSize(1); paramBoxFunc->Draw(); delete [] line; } core->output->setStandardTBOutputText(canvas, ""); core->output->drawAndSaveCanvas(canvas); core->output->tOutFile->cd(Residuals::name.c_str()); h_residualsClone->Write(); delete func; delete funcBox; delete funcGaus; delete h_residualsClone; //delete fitFunc; return; } double fitFunc(double* v, double* par) { //Two error functions to fit a box double x=v[0]; double sigma=par[0]; double area=par[1]; double x0=par[2]; double width=par[3]; double angCoeff=par[4]; double sfita = 0.5*(1+TMath::Erf((x+width-x0)/(TMath::Sqrt(2.)*sigma))); double sfitb = 0.5*(1+TMath::Erf((x-width-x0)/(TMath::Sqrt(2.)*sigma))); return 0.5*(1 + angCoeff*(x-x0))*area*(sfita -sfitb)/(width); } double fitFuncBox(double* v, double* par) { //Two error functions to fit a box double x=v[0]; double sigma=par[0]; double area=par[1]; double width=par[2]; double sfita = 0.5*(1+TMath::Erf((x+width)/(TMath::Sqrt(2.)*sigma))); double sfitb = 0.5*(1+TMath::Erf((x-width)/(TMath::Sqrt(2.)*sigma))); return 0.5*(1 + area*(sfita -sfitb)/(width)); } double multiG(double* v, double* par) { double x = v[0]; double coreMean = par[0]; double outMean = par[1]; double coreWidth = par[2]; double outWidth = par[3]; double coreFrac = par[4]; double norm = par[5]; // G = norm * [coreFrac * C + (1 - coreFrac) * O ] // C: core gaussian double coreArg = (x-coreMean); coreArg /= coreWidth; coreArg *= coreArg; coreWidth /= 2.0; double C = TMath::Exp(-coreArg); C /= (TMath::Sqrt(2.0*TMath::Pi())*coreWidth); // O: outlier gaussian double outArg = (x-outMean); outArg /= outWidth; outArg *= outArg; outWidth /= 2.0; double O = TMath::Exp(-outArg); O /= (TMath::Sqrt(2.0*TMath::Pi())*outWidth); // G = norm *[ coreFrac * C + (1 - coreFrac) * O ] double mG = norm*(coreFrac * C + (1. - coreFrac)* O ); return mG; }