diff --git a/cfg/CoincidenceArrayCo60.cfg b/cfg/CoincidenceArrayCo60.cfg new file mode 100644 index 00000000..0ed63d49 --- /dev/null +++ b/cfg/CoincidenceArrayCo60.cfg @@ -0,0 +1,29 @@ + +inputDirData /data/TOFPET2/reco/ +fileBaseNameData run +inputDirStep1 /data/ALIO/plots/ +inputDirStep2 /data/ALIO/plots/ +fileBaseName1 moduleCharacterization_step +fileBaseName2 run +runs 2506-2517 + + + +outFileName /data/ALIO/plots/coincidenceArrayCo60.root +plotDir /var/www/html/ModuleCharacterization/ALIO/Co60/FBK_thinQuartz/ + + + +tResMin 0 +tResMax 250 +tResMode 2 + +Vov 3 5 7 +energyBins 800 600 400 +energyMin 0 +energyMax 40 + +refBar 7 +refVov 7 +refTh 20 + diff --git a/cfg/moduleCharacterization.cfg b/cfg/moduleCharacterization.cfg index 6fd6dacb..e17603e5 100644 --- a/cfg/moduleCharacterization.cfg +++ b/cfg/moduleCharacterization.cfg @@ -1,19 +1,18 @@ inputDir /data/TOFPET2/reco/ fileBaseName run -runs 2436 +runs 2546 maxEntries -1 -step1FileName /home/cmsdaq/martina/Lab5015Analysis/plots/moduleCharacterization_step1_run2436.root +step1FileName /data/ALIO/plots/moduleCharacterization_step1_run2546.root usePedestals 0 sourceName Na22 -outFileNameStep1 /home/cmsdaq/martina/Lab5015Analysis/plots/moduleCharacterization_step1_run2436.root -outFileNameStep2 /home/cmsdaq/martina/Lab5015Analysis/plots/moduleCharacterization_step2_run2436.root -#plotDir /var/www/html/ModuleCharacterization/FBK_thinQuartz/ -#plotDir /var/www/html/ModuleCharacterization/FBK_siliconeResin/ -plotDir /var/www/html/ModuleCharacterization/HPK_HDR2/ +outFileNameStep1 /data/ALIO/plots/moduleCharacterization_step1_run2546.root +outFileNameStep2 /data/ALIO/plots/moduleCharacterization_step2_run2546.root +plotDir /var/www/html/ModuleCharacterization/alio_guglielmi/corrected/FBK_siliconeResin/ + @@ -23,10 +22,11 @@ tResMax 250 tResMode 2 Vov 3 5 7 -energyBins 800 600 400 +energyBinsNa22 800 600 400 +energyBinsCo60 1000 750 500 energyMin 0 -energyMax 40 - +energyMaxNa22 40 +energyMaxCo60 50 refBar 7 refVov 7 refTh 20 diff --git a/cfg/moduleCharacterization_step3.cfg b/cfg/moduleCharacterization_step3.cfg index 37a678b1..b77c7e17 100644 --- a/cfg/moduleCharacterization_step3.cfg +++ b/cfg/moduleCharacterization_step3.cfg @@ -1,5 +1,6 @@ -inputDir /home/cmsdaq/martina/Lab5015Analysis/plots/ +inputDir1 /home/cmsdaq/martina/Lab5015Analysis/plots/ +inputDir /home/cmsdaq/alio_guglielmi_corr/Lab5015Analysis/plots/ fileBaseName1 moduleCharacterization_step fileBaseName2 run #runs 2428-2453 @@ -8,12 +9,12 @@ runs 2375-2401 -#outFileNameStep3 /home/cmsdaq/martina/Lab5015Analysis/plots/moduleCharacterization_FBK_thinQuartz_step3.root -#plotDirStep3 /var/www/html/ModuleCharacterization/FBK_thinQuartz/ -#outFileNameStep3 /home/cmsdaq/martina/Lab5015Analysis/plots/moduleCharacterization_FBK_siliconeResin_step3.root -#plotDirStep3 /var/www/html/ModuleCharacterization/FBK_siliconeResin/ -outFileNameStep3 /home/cmsdaq/martina/Lab5015Analysis/plots/moduleCharacterization_HPK_HDR2_step3.root -plotDirStep3 /var/www/html/ModuleCharacterization/HPK_HDR2/ +#outFileNameStep3 /home/cmsdaq/alio_guglielmi_corr/Lab5015Analysis/plots/moduleCharacterization_FBK_thinQuartz_step3.root +#plotDirStep3 /var/www/html/ModuleCharacterization/alio_guglielmi/FBK_thinQuartz/ +#outFileNameStep3 /home/cmsdaq/alio_guglielmi_corr/Lab5015Analysis/plots/moduleCharacterization_FBK_siliconeResin_step3.root +#plotDirStep3 /var/www/html/ModuleCharacterization/alio_guglielmi/FBK_siliconeResin/ +outFileNameStep3 /home/cmsdaq/alio_guglielmi_corr/Lab5015Analysis/plots/moduleCharacterization_HPK_HDR2_step3.root +plotDirStep3 /var/www/html/ModuleCharacterization/alio_guglielmi/HPK_HDR2/ diff --git a/main/coincidenceArrayCo60.cpp b/main/coincidenceArrayCo60.cpp new file mode 100644 index 00000000..22da7a31 --- /dev/null +++ b/main/coincidenceArrayCo60.cpp @@ -0,0 +1,1107 @@ +#include "interface/AnalysisUtils.h" +#include "interface/Na22SpectrumAnalyzer.h" +#include "interface/FitUtils.h" +#include "interface/SetTDRStyle.h" +#include "CfgManager/interface/CfgManager.h" +#include "CfgManager/interface/CfgManagerT.h" + +#include +#include +#include +#include +#include +#include +#include + +#include "TFile.h" +#include "TChain.h" +#include "TH1F.h" +#include "TProfile.h" +#include "TProfile2D.h" +#include "TGraphErrors.h" +#include "TF1.h" +#include "TCanvas.h" +#include "TLatex.h" +#include "TLine.h" +#include "TRandom3.h" +#include "TLegend.h" +#include "TSpectrum.h" +#include "TMultiGraph.h" + + + +int main(int argc, char** argv) +{ + setTDRStyle(); + + + if( argc < 2 ) + { + std::cout << ">>> coincideceArrayCo60::usage: " << argv[0] << " configFile.cfg" << std::endl; + return -1; + } + + + //--- parse the config file + CfgManager opts; + opts.ParseConfigFile(argv[1]); + + int debugMode = 0; + if( argc > 2 ) debugMode = atoi(argv[2]); + + + typedef std::numeric_limits dbl; + std::cout.precision(dbl::max_digits10); + + //--- get parameters + std::string plotDir = opts.GetOpt("Output.plotDir"); + //system(Form("rm -r %s/coincidenceResults", plotDir.c_str())); + system(Form("mkdir -p %s",plotDir.c_str())); + system(Form("mkdir -p %s/coincidenceResults/",plotDir.c_str())); + //system(Form("mkdir -p %s/coincidenceResults/tot/",plotDir.c_str())); + system(Form("mkdir -p %s/coincidenceResults/energy/",plotDir.c_str())); + system(Form("mkdir -p %s/coincidenceResults/energy/runs",plotDir.c_str())); + system(Form("mkdir -p %s/coincidenceResults/energy/total",plotDir.c_str())); + system(Form("mkdir -p %s/coincidenceResults/timeResolution/",plotDir.c_str())); + system(Form("mkdir -p %s/coincidenceResults/timeResolution/runs",plotDir.c_str())); + system(Form("mkdir -p %s/coincidenceResults/timeResolution/total",plotDir.c_str())); + system(Form("mkdir -p %s/coincidenceResults/timeResolution/ratioEnergyCorrection",plotDir.c_str())); + system(Form("mkdir -p %s/coincidenceResults/timeResolution/timeResolutionComparison",plotDir.c_str())); + + + //copia il file .php che si trova in una cartella fuori plotDir in plotDir definita nel config file + system(Form("cp %s/../index.php %s",plotDir.c_str(),plotDir.c_str())); + system(Form("cp %s/../index.php %s/coincidenceResults/",plotDir.c_str(),plotDir.c_str())); + //system(Form("cp %s/../index.php %s/coincidenceResults/tot/",plotDir.c_str(),plotDir.c_str())); + system(Form("cp %s/../index.php %s/coincidenceResults/energy/",plotDir.c_str(),plotDir.c_str())); + system(Form("cp %s/../index.php %s/coincidenceResults/energy/runs",plotDir.c_str(),plotDir.c_str())); + system(Form("cp %s/../index.php %s/coincidenceResults/energy/total",plotDir.c_str(),plotDir.c_str())); + system(Form("cp %s/../index.php %s/coincidenceResults/timeResolution/",plotDir.c_str(),plotDir.c_str())); + system(Form("cp %s/../index.php %s/coincidenceResults/timeResolution/runs",plotDir.c_str(),plotDir.c_str())); + system(Form("cp %s/../index.php %s/coincidenceResults/timeResolution/total",plotDir.c_str(),plotDir.c_str())); + system(Form("cp %s/../index.php %s/coincidenceResults/timeResolution/ratioEnergyCorrection",plotDir.c_str(),plotDir.c_str())); + system(Form("cp %s/../index.php %s/coincidenceResults/timeResolution/timeResolutionComparison",plotDir.c_str(),plotDir.c_str())); + + //--- open files and make the tree chain + //--- Prima di eseguire cambia il numero del run nel config file + std::string inputDirStep1 = opts.GetOpt("Input.inputDirStep1"); + std::string inputDirStep2 = opts.GetOpt("Input.inputDirStep2"); + std::string fileBaseName1 = opts.GetOpt("Input.fileBaseName1"); + std::string fileBaseName2 = opts.GetOpt("Input.fileBaseName2"); + std::string inputDirData = opts.GetOpt("Input.inputDirData"); + std::string fileBaseNameData = opts.GetOpt("Input.fileBaseNameData"); + std::string runs = opts.GetOpt("Input.runs"); + + + //Define variables + std::vector listStep1; + std::vector listStep2; + std::vector listData; + std::vector listRun; + + std::stringstream ss(runs); + std::string token; + while( std::getline(ss,token,',') ) + { + std::stringstream ss2(token); + std::string token2; + int runMin = -1; + int runMax = -1; + while( std::getline(ss2,token2,'-') ) + { + if( runMin != -1 && runMax == -1 ) runMax = atoi(token2.c_str()); + if( runMin == -1 ) runMin = atoi(token2.c_str()); + } + if( runMax == -1 ) runMax = runMin; + + for(int run = runMin; run <= runMax; ++run) + { + std::string fileNameStep1; + std::string fileNameStep2; + std::string fileNameData; + fileNameStep1 = Form("%s%s1_%s%04d.root",inputDirStep1.c_str(),fileBaseName1.c_str(),fileBaseName2.c_str(),run); + fileNameStep2 = Form("%s%s2_%s%04d.root",inputDirStep2.c_str(),fileBaseName1.c_str(),fileBaseName2.c_str(),run); + fileNameData = Form("%s/%s%04d_t.root",inputDirData.c_str(),fileBaseNameData.c_str(),run); + std::cout << ">>> Adding file " << fileNameStep1 << std::endl; + std::cout << ">>> Adding file " << fileNameStep2 << std::endl; + std::cout << ">>> Adding file " << fileNameData << std::endl; + + listStep1.push_back(fileNameStep1); + listStep2.push_back(fileNameStep2); + listData.push_back(fileNameData); + listRun.push_back(run); + } + } + + //--- define output root file + std::string outFileName = opts.GetOpt("Output.outFileName"); + TFile* outFile = TFile::Open(outFileName.c_str(),"RECREATE"); + + + std::map trees; + std::map dataTrees; + std::map trees2; + std::vector stepLabels; + std::map< std::string, std::vector > energyRanges; + std::map< std::string, std::map> tProfileFitFunc; + std::map< std::string, std::map> energyRatioFitFunc; + std::map< std::string, std::map> c_timeRes_summary; + std::map< std::string, std::map> l_timeRes_summary; + std::map< std::string, std::map>> gr_timeRes_summary; + + + + //open Step2 root file and data root file + int id = 0; + + for(auto file: listStep2){ + std::string label; + std::string labelData; + std::string label2; + TFile* inFile = TFile::Open(file.c_str(),"READ"); + + TList* list = inFile -> GetListOfKeys(); + TIter next(list); + TObject* object = 0; + while( (object = next()) ){ + std::string name(object->GetName()); + std::vector tokens = GetTokens(name,'_'); + std::size_t found; + std::size_t found2; + found = name.find("data_"); + found2 = name.find("dataRes_"); + if( found!=std::string::npos ){ + label = Form("%s_%s_%s_%04d", tokens[1].c_str(), tokens[2].c_str(), tokens[3].c_str(), listRun[id]); + labelData = Form("%s_%s_%04d", tokens[2].c_str(), tokens[3].c_str(), listRun[id]); + trees[label] = (TTree*)( inFile->Get(name.c_str()) ); + } + if( found2!=std::string::npos ){ + label2 = Form("%s_%s_%s_%s_%04d", tokens[1].c_str(), tokens[2].c_str(), tokens[3].c_str(), tokens[4].c_str(), listRun[id]); + trees2[label2] = (TTree*)( inFile->Get(name.c_str()) ); + } + } + dataTrees[labelData] = new TChain("data","data"); + dataTrees[labelData] -> Add(listData[id].c_str()); + id ++; + } + + + + + //taking 7-8-9 energyBin from filled Trees + + + for (auto mapIt : trees){ + + float energyBinVal[3]; + float theIndex = -1; + mapIt.second -> SetBranchStatus("*",0); + mapIt.second -> SetBranchStatus("energyBinValCo60", 1); mapIt.second -> SetBranchAddress("energyBinValCo60", energyBinVal); + mapIt.second -> SetBranchStatus("indexID", 1); mapIt.second -> SetBranchAddress("indexID", &theIndex); + int nEntries = mapIt.second->GetEntries(); + for(int entry = 0; entry < nEntries; ++entry){ + mapIt.second -> GetEntry(entry); + + std::vector vec; + if (theIndex > 0){ + for ( int i = 0; i < 3; i++){ + vec.push_back(energyBinVal[i]); + } + std::sort(vec.begin(), vec.end()); + energyRanges[mapIt.first] = vec; + } + + } + + } + + for ( auto mapIt : trees2){ + int iBin; + float tRes; + int run; + mapIt.second -> SetBranchStatus("*",0); + mapIt.second -> SetBranchStatus("timeResolution", 1); mapIt.second -> SetBranchAddress("timeResolution", &tRes); + int nEntries = mapIt.second->GetEntries(); + std::vector tokens = GetTokens(mapIt.first,'_'); + std::string label_summary = Form( "%s_%s_%s",tokens[0].c_str(), tokens[1].c_str(), tokens[2].c_str()); + std::vector tokens2 = GetTokens(tokens[3].c_str(),'n'); + iBin = int( atof(tokens2[2].c_str())); + run = int( atof(tokens[4].c_str())); + for(int entry = 0; entry < nEntries; ++entry){ + mapIt.second -> GetEntry(entry); + + if ( iBin == 8 || iBin == 9){ + if (c_timeRes_summary[label_summary][iBin] == NULL){ + c_timeRes_summary[label_summary][iBin] = new TCanvas ( Form("c_time_resolution_comparison%s_energyBin%d", label_summary.c_str(), iBin), Form("c_time_resolution_comparison%s_energyBin%d", label_summary.c_str(), iBin)); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> Draw(); + gPad -> SetGridy(); + l_timeRes_summary[label_summary][iBin] = new TLegend (0.70,0.40,0.95,0.95); + l_timeRes_summary[label_summary][iBin] -> SetBorderSize(0); + } + + TGraphErrors* graph = new TGraphErrors(); + graph -> SetPoint (0, gr_timeRes_summary[label_summary][iBin].size()+1, tRes); + graph -> SetPointError(0, 0, 0); + graph -> SetMarkerColor(1 + gr_timeRes_summary[label_summary][iBin].size()); + graph -> SetMarkerStyle(25); + gr_timeRes_summary[label_summary][iBin].push_back((TGraphErrors*) graph); + l_timeRes_summary[label_summary][iBin] -> AddEntry ( graph, Form("run %4d", run), "p"); + } + } + } + + + //bars of interest + std::map > barMap; + for( auto index :energyRanges){ + std::vector tokens = GetTokens(index.first,'_'); + std::string label = Form("%s_%s_%s", tokens[1].c_str(), tokens[2].c_str(), tokens[3].c_str()); + std::vector tokens2 = GetTokens(tokens[0],'r'); + std::vector tokens3 = GetTokens(tokens2[1],'L'); + int iBar = int(atof(tokens3[0].c_str())); + barMap[label].push_back(iBar); + } + + //taking TProfiles and energyRatio for 8-9 energyBin and bars of interest + id = 0; + std::map < std::string, TCanvas*> c_step2_TProfile_run; + for(auto file: listStep2){ + std::string label; + std::string labelTProfile; + std::string labelEnergyRatio; + TFile* inFile = TFile::Open(file.c_str(),"READ"); + + for ( auto mapIt : energyRanges){ + label = mapIt.first; + std::vector tokens = GetTokens(mapIt.first,'_'); + if ( atof(tokens[3].c_str()) != listRun[id]) continue; + for( int iBin = 8; iBin < 10; iBin++){ + std::string labelEnFile = Form("h1_energyRatio_%s_%s_%s_energyBin%d", tokens[0].c_str(), tokens[1].c_str(), tokens[2].c_str(), iBin); + labelEnergyRatio = Form("h1_energyRatio_%s_%s_%s_energyBin%d_%04d", tokens[0].c_str(), tokens[1].c_str(), tokens[2].c_str(), iBin, listRun[id]); + TCanvas * canvas = new TCanvas ( Form("canvas_%s",labelEnergyRatio.c_str() ), Form("canvas_%s",labelEnergyRatio.c_str() )); + canvas -> cd(); + + TH1F * histo = (TH1F*) ( inFile -> Get( Form("%s", labelEnFile.c_str()))); + if( !histo) continue; + energyRatioFitFunc[label][iBin] = new TF1(Form("energyRatioFitFunc%s",Form("bin%d_%s", iBin, label.c_str())),"gaus",histo->GetMean()-2.*histo->GetRMS(),histo->GetMean()+2.*histo->GetRMS()); + histo -> Fit(energyRatioFitFunc[label][iBin],"QNRS"); + histo -> Fit(energyRatioFitFunc[label][iBin],"QS+","",histo->GetMean()-2.*histo->GetRMS(),histo->GetMean()+2.*histo->GetRMS()); + outFile -> cd(); + histo -> Write(); + + + + labelTProfile = Form("p1_deltaT_vs_energyRatio_%s_%s_%s_energyBin%d", tokens[0].c_str(), tokens[1].c_str(), tokens[2].c_str(), iBin); + std::string labelCanvas = Form( "%s_%s_%s_energyBin%d_%04d",tokens[0].c_str(), tokens[1].c_str(), tokens[2].c_str(), iBin, listRun[id]); + + TProfile* prof = (TProfile*)( inFile->Get(Form("%s",labelTProfile.c_str()))); + if( !prof) continue; + c_step2_TProfile_run[labelCanvas] = new TCanvas(Form("c_deltaT_vs_energyRatio_%s", labelCanvas.c_str()), Form("c_deltaT_vs_energyRatio_%s", labelCanvas.c_str())); + + tProfileFitFunc[label][iBin] = new TF1(Form("tProfileFitFunc%s",Form("bin%d_%s", iBin, label.c_str())),"pol1",-0.88,1.12); + prof -> Fit(tProfileFitFunc[label][iBin],"QRS+"); + + + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> Draw(); + gPad -> SetGridy(); + prof -> SetTitle(Form(";energy_{right} / energy_{left};#Deltat [ps]")); + prof -> Draw(""); + tProfileFitFunc[label][iBin] -> SetLineColor(kRed); + tProfileFitFunc[label][iBin] -> SetLineWidth(2); + tProfileFitFunc[label][iBin] -> Draw("same"); + outFile -> cd(); + prof -> Write(); + + + } + } + id++; + } + + + + + + int chsL[16]; + int chsR[16]; + chsL[0] = 249; + chsR[0] = 4; + chsL[1] = 242; + chsR[1] = 6; + chsL[2] = 251; + chsR[2] = 15; + chsR[3] = 8; + chsL[4] = 250; + chsR[4] = 13; + chsL[5] = 254; + chsR[5] = 2; + chsL[6] = 201; + chsR[6] = 11; + chsL[7] = 230; + chsR[7] = 27; + chsL[8] = 226; + chsR[8] = 32; + chsL[9] = 225; + chsR[9] = 31; + chsL[10] = 228; + chsR[10] = 30; + chsL[11] = 227; + chsR[11] = 29; + chsL[12] = 229; + chsR[12] = 28; + chsL[13] = 231; + chsR[13] = 26; + chsL[14] = 233; + chsR[14] = 24; + chsL[15] = 232; + chsR[15] = 25; + + std::map < std::string, TCanvas*> c_energy_bar_run; + std::map < std::string, TCanvas*> c_energy_bar_tot; + std::map < std::string, TCanvas*> c_energy_1173_bar_run; + std::map < std::string, TCanvas*> c_energy_1173_bar_tot; + std::map < std::string, TCanvas*> c_energy_1332_bar_run; + std::map < std::string, TCanvas*> c_energy_1332_bar_tot; + std::map < std::string, TCanvas*> c_energy_coincidence_run; + std::map < std::string, TCanvas*> c_energy_coincidence_tot; + std::map < std::string, TCanvas*> c_timeRes_coincidence_run; + std::map < std::string, TCanvas*> c_timeRes_coincidence_tot; + std::map < std::string, std::map > c_TProfile_energyRatio_corr; + + std::map < std::string, TH1F*> h_energy_bar_run; + std::map < std::string, TH1F*> h_energy_bar_tot; + std::map < std::string, TH1F*> h_energy_1173_bar_run; + std::map < std::string, TH1F*> h_energy_1173_bar_tot; + std::map < std::string, TH1F*> h_energy_1332_bar_run; + std::map < std::string, TH1F*> h_energy_1332_bar_tot; + std::map < std::string, TH1F*> h_energy_coincidence_run; + std::map < std::string, TH1F*> h_energy_coincidence_tot; + std::map < std::string, TH1F*> h_timeRes_coincidence_run; + std::map < std::string, TH1F*> h_timeRes_coincidence_tot; + std::map < std::string, TH1F*> h_timeRes_corr_coincidence_run; + std::map < std::string, TH1F*> h_timeRes_corr_coincidence_tot; + std::map < std::string, TH1F*> h_timeRes_corr_data_coincidence_run; + std::map < std::string, TH1F*> h_timeRes_corr_data_coincidence_tot; + std::map < std::string, std::map > p1_TProfile_energyRatio_corr; + std::map < std::string, std::map > h1_energyRatio; + + for (auto mapIt : dataTrees){ + for( int iBar : barMap[mapIt.first]){ + + std::string label_c = Form( "bar%02dL-R_%s", iBar, mapIt.first.c_str()); + std::vector tokens = GetTokens(mapIt.first,'_'); + std::string labelTot_c = Form( "bar%02dL-R_%s_%s", iBar, tokens[0].c_str(), tokens[1].c_str()); + + + + //all spectrum histo and canvan + if(c_energy_bar_run[label_c]==NULL){ + c_energy_bar_run[label_c] = new TCanvas(Form("c_energy_%s",label_c.c_str()), Form("c_energy_%s",label_c.c_str())); + c_energy_bar_run[label_c] -> cd(); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> SetTitle(";entries;energy [a.u.]"); + hPad -> Draw(); + gPad -> SetGridy(); + h_energy_bar_run[label_c] = new TH1F ( Form("h_energy_%s",label_c.c_str()), Form("h_energy_%s",label_c.c_str()), 1000, 0, 50); + } + + if(c_energy_bar_tot[labelTot_c]==NULL){ + c_energy_bar_tot[labelTot_c] = new TCanvas(Form("c_energy_%s",labelTot_c.c_str()), Form("c_energy_%s",labelTot_c.c_str())); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> SetTitle(";entries;energy [a.u.]"); + hPad -> Draw(); + gPad -> SetGridy(); + h_energy_bar_tot[labelTot_c] = new TH1F ( Form("h_energy_%s",labelTot_c.c_str()), Form("h_energy_%s",labelTot_c.c_str()), 1000, 0, 50); + } + //1173 peak histo and canvan + if(c_energy_1173_bar_run[label_c]==NULL){ + c_energy_1173_bar_run[label_c] = new TCanvas(Form("c_1173_peak_%s",label_c.c_str()), Form("c_1173_peak_%s",label_c.c_str())); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> SetTitle(";entries;energy [a.u.]"); + hPad -> Draw(); + gPad -> SetGridy(); + h_energy_1173_bar_run[label_c] = new TH1F ( Form("h_1173_peak_%s",label_c.c_str()), Form("h_1173_peak_%s",label_c.c_str()), 1000, 0, 50); + } + + if(c_energy_1173_bar_tot[labelTot_c]==NULL){ + c_energy_1173_bar_tot[labelTot_c] = new TCanvas(Form("c_1173_peak_%s",labelTot_c.c_str()), Form("c_1173_peak_%s",labelTot_c.c_str())); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> SetTitle(";entries;energy [a.u.]"); + hPad -> Draw(); + gPad -> SetGridy(); + h_energy_1173_bar_tot[labelTot_c] = new TH1F ( Form("h_1173_peak_%s",labelTot_c.c_str()), Form("h_1173_peak_%s",labelTot_c.c_str()), 1000, 0, 50); + } + //1332 peak histo and canvan + if(c_energy_1332_bar_run[label_c]==NULL){ + c_energy_1332_bar_run[label_c] = new TCanvas(Form("c_1332_peak_%s",label_c.c_str()), Form("c_1332_peak_%s",label_c.c_str())); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> SetTitle(";entries;energy [a.u.]"); + hPad -> Draw(); + gPad -> SetGridy(); + h_energy_1332_bar_run[label_c] = new TH1F ( Form("h_1332_peak_%s",label_c.c_str()), Form("h_1332_peak_%s",label_c.c_str()), 1000, 0, 50); + } + + if(c_energy_1332_bar_tot[labelTot_c]==NULL){ + c_energy_1332_bar_tot[labelTot_c] = new TCanvas(Form("c_1332_peak_%s",labelTot_c.c_str()), Form("c_1332_peak_%s",labelTot_c.c_str())); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> SetTitle(";entries;energy [a.u.]"); + hPad -> Draw(); + gPad -> SetGridy(); + h_energy_1332_bar_tot[labelTot_c] = new TH1F ( Form("h_1332_peak_%s",labelTot_c.c_str()), Form("h_1332_peak_%s",labelTot_c.c_str()), 1000, 0, 50); + } + //peaks -> iBar coincidence with the others (energy and time resolution) + for (int coinBar : barMap[mapIt.first]){ + if(coinBar != iBar){ + std::string label_coin = Form( "1173_bar%02dL-R_vs_1332_bar%02dL-R_%s", iBar, coinBar, mapIt.first.c_str()); + std::string labelTot_coin = Form( "1173_bar%02dL-R_vs_1332_bar%02dL-R_%s_%s", iBar, coinBar, tokens[0].c_str(), tokens[1].c_str()); + std::string label_coin2 = Form( "1332_bar%02dL-R_vs_1173_bar%02dL-R_%s", iBar, coinBar, mapIt.first.c_str()); + std::string labelTot_coin2 = Form( "1332_bar%02dL-R_vs_1173_bar%02dL-R_%s_%s", iBar, coinBar, tokens[0].c_str(), tokens[1].c_str()); + + if(c_energy_coincidence_run[label_coin]==NULL){ + c_energy_coincidence_run[label_coin] = new TCanvas(Form("c_energy_coincidence_%s",label_coin.c_str()), Form("c_energy_coincidence_%s",label_coin.c_str())); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> SetTitle(";entries;energy [a.u.]"); + hPad -> Draw(); + gPad -> SetGridy(); + h_energy_coincidence_run[label_coin] = new TH1F ( Form("h_energy_coincidence_%s",label_coin.c_str()), Form("h_energy_coincidence_%s",label_coin.c_str()), 1000, 0, 50); + } + + if(c_energy_coincidence_tot[labelTot_coin]==NULL){ + c_energy_coincidence_tot[labelTot_coin] = new TCanvas(Form("c_energy_coincidence_%s",labelTot_coin.c_str()), Form("c_energy_coincidence_%s",labelTot_coin.c_str())); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> SetTitle(";entries;energy [a.u.]"); + hPad -> Draw(); + gPad -> SetGridy(); + h_energy_coincidence_tot[labelTot_coin] = new TH1F ( Form("h_energy_coincidence_%s",labelTot_coin.c_str()), Form("h_energy_coincidence_%s",labelTot_coin.c_str()), 1000, 0, 50); + } + + if(c_energy_coincidence_run[label_coin2]==NULL){ + c_energy_coincidence_run[label_coin2] = new TCanvas(Form("c_energy_coincidence_%s",label_coin2.c_str()), Form("c_energy_coincidence_%s",label_coin2.c_str())); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> SetTitle(";entries;energy [a.u.]"); + hPad -> Draw(); + gPad -> SetGridy(); + h_energy_coincidence_run[label_coin2] = new TH1F ( Form("h_energy_coincidence_%s",label_coin2.c_str()), Form("h_energy_coincidence_%s",label_coin2.c_str()), 1000, 0, 50); + } + + if(c_energy_coincidence_tot[labelTot_coin2]==NULL){ + c_energy_coincidence_tot[labelTot_coin2] = new TCanvas(Form("c_energy_coincidence_%s",labelTot_coin2.c_str()), Form("c_energy_coincidence_%s",labelTot_coin2.c_str())); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> SetTitle(";entries;energy [a.u.]"); + hPad -> Draw(); + gPad -> SetGridy(); + h_energy_coincidence_tot[labelTot_coin2] = new TH1F ( Form("h_energy_coincidence_%s",labelTot_coin2.c_str()), Form("h_energy_coincidence_%s",labelTot_coin2.c_str()), 1000, 0, 50); + } + if(c_timeRes_coincidence_run[label_coin]==NULL){ + c_timeRes_coincidence_run[label_coin] = new TCanvas(Form("c_timeRes_coincidence_%s",label_coin.c_str()), Form("c_timeRes_coincidence_%s",label_coin.c_str())); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> SetTitle("energy-corrected #Deltat [ps];entries"); + hPad -> Draw(); + gPad -> SetGridy(); + h_timeRes_coincidence_run[label_coin] = new TH1F ( Form("h_timeRes_coincidence_%s",label_coin.c_str()), Form("h_timeRes_coincidence_%s",label_coin.c_str()), 200, -4000, 4000); + h_timeRes_corr_coincidence_run[label_coin] = new TH1F ( Form("h_timeRes_corr_coincidence_%s",label_coin.c_str()), Form("h_timeRes_corr_coincidence_%s",label_coin.c_str()), 200, -4000, 4000); + h_timeRes_corr_data_coincidence_run[label_coin] = new TH1F ( Form("h_timeRes_corr_data_coincidence_%s",label_coin.c_str()), Form("h_timeRes_corr_data_coincidence_%s",label_coin.c_str()), 200, -4000, 4000); + } + + if(c_timeRes_coincidence_tot[labelTot_coin]==NULL){ + c_timeRes_coincidence_tot[labelTot_coin] = new TCanvas(Form("c_timeRes_coincidence_%s",labelTot_coin.c_str()), Form("c_timeRes_coincidence_%s",labelTot_coin.c_str())); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> SetTitle(";energy-corrected #Deltat [ps];entries"); + hPad -> Draw(); + gPad -> SetGridy(); + h_timeRes_coincidence_tot[labelTot_coin] = new TH1F ( Form("h_timeRes_coincidence_%s",labelTot_coin.c_str()), Form("h_timeRes_coincidence_%s",labelTot_coin.c_str()), 200, -4000, 4000); + h_timeRes_corr_coincidence_tot[labelTot_coin] = new TH1F ( Form("h_timeRes_corr_coincidence_%s",labelTot_coin.c_str()), Form("h_timeRes_corr_coincidence_%s",labelTot_coin.c_str()), 200, -4000, 4000); + h_timeRes_corr_data_coincidence_tot[labelTot_coin] = new TH1F ( Form("h_timeRes_corr_data_coincidence_%s",labelTot_coin.c_str()), Form("h_timeRes_corr_data_coincidence_%s",labelTot_coin.c_str()), 200, -4000, 4000); + } + + if(c_timeRes_coincidence_run[label_coin2]==NULL){ + c_timeRes_coincidence_run[label_coin2] = new TCanvas(Form("c_timeRes_coincidence_%s",label_coin2.c_str()), Form("c_timeRes_coincidence_%s",label_coin2.c_str())); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> SetTitle("energy-corrected #Deltat [ps];entries"); + hPad -> Draw(); + gPad -> SetGridy(); + h_timeRes_coincidence_run[label_coin2] = new TH1F ( Form("h_timeRes_coincidence_%s",label_coin2.c_str()), Form("h_timeRes_coincidence_%s",label_coin2.c_str()), 200, -4000, 4000); + h_timeRes_corr_coincidence_run[label_coin2] = new TH1F ( Form("h_timeRes_corr_coincidence_%s",label_coin2.c_str()), Form("h_timeRes_corr_coincidence_%s",label_coin2.c_str()), 200, -4000, 4000); + h_timeRes_corr_data_coincidence_run[label_coin2] = new TH1F ( Form("h_timeRes_corr_data_coincidence_%s",label_coin2.c_str()), Form("h_timeRes_corr_data_coincidence_%s",label_coin2.c_str()), 200, -4000, 4000); + } + + if(c_timeRes_coincidence_tot[labelTot_coin2]==NULL){ + c_timeRes_coincidence_tot[labelTot_coin2] = new TCanvas(Form("c_timeRes_coincidence_%s",labelTot_coin2.c_str()), Form("c_timeRes_coincidence_%s",labelTot_coin2.c_str())); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> SetTitle(";energy-corrected #Deltat [ps];entries"); + hPad -> Draw(); + gPad -> SetGridy(); + h_timeRes_coincidence_tot[labelTot_coin2] = new TH1F ( Form("h_timeRes_coincidence_%s",labelTot_coin2.c_str()), Form("h_timeRes_coincidence_%s",labelTot_coin2.c_str()), 200, -4000, 4000); + h_timeRes_corr_coincidence_tot[labelTot_coin2] = new TH1F ( Form("h_timeRes_corr_coincidence_%s",labelTot_coin2.c_str()), Form("h_timeRes_corr_coincidence_%s",labelTot_coin2.c_str()), 200, -4000, 4000); + h_timeRes_corr_data_coincidence_tot[labelTot_coin2] = new TH1F ( Form("h_timeRes_corr_data_coincidence_%s",labelTot_coin2.c_str()), Form("h_timeRes_corr_data_coincidence_%s",labelTot_coin2.c_str()), 200, -4000, 4000); + } + } + } + for(int iBin = 8; iBin < 10; iBin ++){ + if ( c_TProfile_energyRatio_corr[labelTot_c][iBin] == NULL){ + c_TProfile_energyRatio_corr[labelTot_c][iBin] = new TCanvas( Form("c_deltaT_vs_energyRatio_%s_energyBin%d", labelTot_c.c_str(), iBin), Form("c_deltaT_vs_energyRatio_%s_energyBin%d", labelTot_c.c_str(),iBin)); + TH1F* hPad = (TH1F*)( gPad->DrawFrame(-1.,0.,64.,25.) ); + hPad -> Draw(); + gPad -> SetGridy(); + p1_TProfile_energyRatio_corr[labelTot_c][iBin] = new TProfile(Form("p1_deltaT_vs_energyRatio_%s_energyBin%d",labelTot_c.c_str(), iBin),Form("p1_deltaT_vs_energyRatio_%s_energyBin%d",labelTot_c.c_str(), iBin),100, 0, 2, -4000,4000); + h1_energyRatio[labelTot_c][iBin] = new TH1F(Form("h1_energyRatio_%s_energyBin%d",labelTot_c.c_str(),iBin),"",1000,0.,5.); + } + } + } + } + + std::map> acceptEvent; + + for (auto mapIt : dataTrees){ + std::cout << mapIt.first << std::endl; + int nEvents = mapIt.second -> GetEntries(); + + float qfine[256]; + float energy[256]; + long long time[256]; + mapIt.second -> SetBranchStatus("*",0); + mapIt.second -> SetBranchStatus("qfine", 1); mapIt.second -> SetBranchAddress("qfine", qfine); + mapIt.second -> SetBranchStatus("energy", 1); mapIt.second -> SetBranchAddress("energy", energy); + mapIt.second -> SetBranchStatus("time", 1); mapIt.second -> SetBranchAddress("time", time); + + for (int entry = 0; entry < nEvents; entry++){ + mapIt.second -> GetEntry(entry); + if (entry%500000 == 0) std::cout << "Analyzing entry " << entry << "/" << nEvents << std::endl; + acceptEvent[mapIt.first][entry] = false; + for( int iBar : barMap[mapIt.first]){ + + std::string label_c = Form( "bar%02dL-R_%s", iBar, mapIt.first.c_str()); + std::string label_en = Form( "bar%02dL-R_%s", iBar, mapIt.first.c_str()); + std::vector tokens = GetTokens(mapIt.first,'_'); + std::string labelTot_c = Form( "bar%02dL-R_%s_%s", iBar, tokens[0].c_str(), tokens[1].c_str()); + + int chL(chsL[iBar]); + int chR(chsR[iBar]); + + float qfineL(qfine[chL]); + float qfineR(qfine[chR]); + + if( qfineL < 13. || qfineR < 13. ) continue; + if( qfineL > 500. || qfineR > 500. ) continue; + + float energyL(energy[chL]); + float energyR(energy[chR]); + long long timeL(time[chL]); + long long timeR(time[chR]); + float avEnergy = 0.5*(energyL+energyR); + //energy spectrum + h_energy_bar_run[label_c] -> Fill( avEnergy); + h_energy_bar_tot[labelTot_c] -> Fill( avEnergy); + + if( avEnergy > energyRanges[label_en][0] && avEnergy < energyRanges[label_en][1] ){ + + //1173 peak + h_energy_1173_bar_run[label_c] -> Fill( avEnergy); + h_energy_1173_bar_tot[labelTot_c] -> Fill( avEnergy); + + //coincidence 1173 vs 1332 (no energy in other bars) + int coin_chL[barMap[mapIt.first].size()]; + int coin_chR[barMap[mapIt.first].size()]; + float coin_qfineL[barMap[mapIt.first].size()]; + float coin_qfineR[barMap[mapIt.first].size()]; + + int counter = 0; + for ( int i = 0; i < barMap[mapIt.first].size(); i++){ + int coinBar = barMap[mapIt.first][i]; + coin_chL[i] = chsL[coinBar]; + coin_chR[i] = chsR[coinBar]; + coin_qfineL[i] = qfine[coin_chL[i]]; + coin_qfineR[i] = qfine[coin_chR[i]]; + if ( coin_qfineL[i] < 13. && coin_qfineR[i] < 13.){ + counter ++; + } + } + if(counter == barMap[mapIt.first].size() -2){ + for ( int i = 0; i < barMap[mapIt.first].size(); i++){ + int coinBar = barMap[mapIt.first][i]; + if (coinBar != iBar && coin_qfineL[i] > 13. && coin_qfineR[i] > 13.){ + float coinEnergy = 0.5 *( energy[coin_chL[i]] + energy[coin_chR[i]]); + if( coinEnergy > energyRanges[label_en][1] && coinEnergy < energyRanges[label_en][2] ){ + std::string label_coin = Form( "1173_bar%02dL-R_vs_1332_bar%02dL-R_%s", iBar, coinBar, mapIt.first.c_str()); + std::string labelTot_coin = Form( "1173_bar%02dL-R_vs_1332_bar%02dL-R_%s_%s", iBar, coinBar, tokens[0].c_str(), tokens[1].c_str()); + h_energy_coincidence_run[label_coin] -> Fill(avEnergy); + h_energy_coincidence_tot[labelTot_coin] -> Fill(avEnergy); + //time resolution + h_timeRes_coincidence_run[label_coin] -> Fill(timeR - timeL); + h_timeRes_coincidence_tot[labelTot_coin] -> Fill(timeR -timeL); + acceptEvent[mapIt.first][entry] = true; + + float energyRatioCorr = tProfileFitFunc[label_en][9]->Eval(energyR/energyL) - tProfileFitFunc[label_en][9]->Eval(energyRatioFitFunc[label_en][9]->GetParameter(1)); + h_timeRes_corr_coincidence_run[label_coin] -> Fill(timeR -timeL -energyRatioCorr); + h_timeRes_corr_coincidence_tot[labelTot_coin] -> Fill(timeR -timeL -energyRatioCorr); + + h1_energyRatio[labelTot_c][8] -> Fill (energyR/energyL); + p1_TProfile_energyRatio_corr[labelTot_c][8] -> Fill( energyR/energyL, timeR - timeL); + } + } + } + } + } + + if( avEnergy > energyRanges[label_en][1] && avEnergy < energyRanges[label_en][2] ){ + + //1332 peak + h_energy_1332_bar_run[label_c] -> Fill( avEnergy); + h_energy_1332_bar_tot[labelTot_c] -> Fill( avEnergy); + + //coincidence 1332 vs 1173 (no energy in other bars) + int coin_chL[barMap[mapIt.first].size()]; + int coin_chR[barMap[mapIt.first].size()]; + float coin_qfineL[barMap[mapIt.first].size()]; + float coin_qfineR[barMap[mapIt.first].size()]; + + int counter = 0; + for ( int i = 0; i < barMap[mapIt.first].size(); i++){ + int coinBar = barMap[mapIt.first][i]; + coin_chL[i] = chsL[coinBar]; + coin_chR[i] = chsR[coinBar]; + coin_qfineL[i] = qfine[coin_chL[i]]; + coin_qfineR[i] = qfine[coin_chR[i]]; + if ( coin_qfineL[i] < 13. && coin_qfineR[i] < 13.){ + counter ++; + } + } + if(counter == barMap[mapIt.first].size() -2){ + for ( int i = 0; i < barMap[mapIt.first].size(); i++){ + int coinBar = barMap[mapIt.first][i]; + if (coinBar != iBar && coin_qfineL[i] > 13. && coin_qfineR[i] > 13.){ + float coinEnergy = 0.5 *( energy[coin_chL[i]] + energy[coin_chR[i]]); + if( coinEnergy > energyRanges[label_en][0] && coinEnergy < energyRanges[label_en][1] ){ + std::string label_coin = Form( "1332_bar%02dL-R_vs_1173_bar%02dL-R_%s", iBar, coinBar, mapIt.first.c_str()); + std::string labelTot_coin = Form( "1332_bar%02dL-R_vs_1173_bar%02dL-R_%s_%s", iBar, coinBar, tokens[0].c_str(), tokens[1].c_str()); + h_energy_coincidence_run[label_coin] -> Fill(avEnergy); + h_energy_coincidence_tot[labelTot_coin] -> Fill(avEnergy); + + //time resolution + h_timeRes_coincidence_run[label_coin] -> Fill(timeR - timeL); + h_timeRes_coincidence_tot[labelTot_coin] -> Fill(timeR -timeL); + acceptEvent[mapIt.first][entry] = true; + + float energyRatioCorr = tProfileFitFunc[label_en][9]->Eval(energyR/energyL) - tProfileFitFunc[label_en][9]->Eval(energyRatioFitFunc[label_en][9]->GetParameter(1)); + h_timeRes_corr_coincidence_run[label_coin] -> Fill(timeR -timeL -energyRatioCorr); + h_timeRes_corr_coincidence_tot[labelTot_coin] -> Fill(timeR -timeL -energyRatioCorr); + h1_energyRatio[labelTot_c][9] -> Fill (energyR/energyL); + p1_TProfile_energyRatio_corr[labelTot_c][9] -> Fill( energyR/energyL, timeR - timeL); + } + } + } + } + } + } + } + } + + + + std::map < std::string, std::map > p1_TProfile_energyRatio_corr_FitFunc; + std::map < std::string, std::map > h1_energyRatio_FitFunc; + + for (auto mapIt :c_TProfile_energyRatio_corr){ + for(int iBin = 8; iBin < 10; iBin++){ + + c_TProfile_energyRatio_corr[mapIt.first][iBin] -> cd(); + TProfile *prof = p1_TProfile_energyRatio_corr[mapIt.first][iBin]; + p1_TProfile_energyRatio_corr_FitFunc[mapIt.first][iBin] = new TF1(Form("p1_TProfile_energyRatio_corr_FitFunc%s",Form("bin%d_%s", iBin, mapIt.first.c_str())),"pol4",0,10); + float xmin = p1_TProfile_energyRatio_corr[mapIt.first][iBin] -> GetMean() - 2.5*p1_TProfile_energyRatio_corr[mapIt.first][iBin] -> GetRMS(); + float xmax = p1_TProfile_energyRatio_corr[mapIt.first][iBin] -> GetMean() + 2.5*p1_TProfile_energyRatio_corr[mapIt.first][iBin] -> GetRMS(); + p1_TProfile_energyRatio_corr_FitFunc[mapIt.first][iBin] -> SetRange(xmin, xmax); + prof -> Fit(p1_TProfile_energyRatio_corr_FitFunc[mapIt.first][iBin],"QRS+"); + prof -> SetTitle(Form(";energy_{right} / energy_{left};#Deltat [ps]")); + prof -> Draw(""); + p1_TProfile_energyRatio_corr_FitFunc[mapIt.first][iBin] -> SetLineColor(kRed); + p1_TProfile_energyRatio_corr_FitFunc[mapIt.first][iBin] -> SetLineWidth(2); + p1_TProfile_energyRatio_corr_FitFunc[mapIt.first][iBin] -> Draw("same"); + + outFile -> cd(); + prof -> Write(); + } + } + for (auto mapIt : h1_energyRatio){ + for(int iBin = 8; iBin < 10; iBin++){ + + TCanvas * canvas = new TCanvas (Form("canvas_energyRatio%s",Form("bin%d_%s", iBin, mapIt.first.c_str())), Form("canvas_energyRatio%s",Form("bin%d_%s", iBin, mapIt.first.c_str()))); + canvas -> cd(); + TH1F* histo = h1_energyRatio[mapIt.first][iBin]; + h1_energyRatio_FitFunc[mapIt.first][iBin] = new TF1(Form("h1_energyRatio_FitFunc%s",Form("bin%d_%s", iBin, mapIt.first.c_str())),"gaus",histo->GetMean()-2.*histo->GetRMS(),histo->GetMean()+2.*histo->GetRMS()); + histo -> Fit(h1_energyRatio_FitFunc[mapIt.first][iBin],"QNRS"); + histo -> Fit(h1_energyRatio_FitFunc[mapIt.first][iBin],"QS+","",histo->GetMean()-2.*histo->GetRMS(),histo->GetMean()+2.*histo->GetRMS()); + + outFile -> cd(); + histo -> Write(); + } + } + + + + std::cout << " Second loop over events " << std::endl; + for (auto mapIt : dataTrees){ + std::cout << mapIt.first << std::endl; + int nEvents = mapIt.second -> GetEntries(); + + float qfine[256]; + float energy[256]; + long long time[256]; + + for (int entry = 0; entry < nEvents; entry++){ + mapIt.second -> GetEntry(entry); + if (entry%500000 == 0) std::cout << "Analyzing entry " << entry << "/" << nEvents << std::endl; + if( !acceptEvent[mapIt.first][entry] ) continue; + for( int iBar : barMap[mapIt.first]){ + + std::string label_c = Form( "bar%02dL-R_%s", iBar, mapIt.first.c_str()); + std::vector tokens = GetTokens(mapIt.first,'_'); + std::string labelTot_c = Form( "bar%02dL-R_%s_%s", iBar, tokens[0].c_str(), tokens[1].c_str()); + + int chL(chsL[iBar]); + int chR(chsR[iBar]); + + float qfineL(qfine[chL]); + float qfineR(qfine[chR]); + + if( qfineL < 13. || qfineR < 13. ) continue; + if( qfineL > 500. || qfineR > 500. ) continue; + + float energyL(energy[chL]); + float energyR(energy[chR]); + long long timeL(time[chL]); + long long timeR(time[chR]); + float avEnergy = 0.5*(energyL+energyR); + + + if( avEnergy > energyRanges[label_c][0] && avEnergy < energyRanges[label_c][1] ){ + + + //coincidence 1173 vs 1332 (no energy in other bars) + int coin_chL[barMap[mapIt.first].size()]; + int coin_chR[barMap[mapIt.first].size()]; + float coin_qfineL[barMap[mapIt.first].size()]; + float coin_qfineR[barMap[mapIt.first].size()]; + + int counter = 0; + for ( int i = 0; i < barMap[mapIt.first].size(); i++){ + int coinBar = barMap[mapIt.first][i]; + coin_chL[i] = chsL[coinBar]; + coin_chR[i] = chsR[coinBar]; + coin_qfineL[i] = qfine[coin_chL[i]]; + coin_qfineR[i] = qfine[coin_chR[i]]; + if ( coin_qfineL[i] < 13. && coin_qfineR[i] < 13.){ + counter ++; + } + } + if(counter == barMap[mapIt.first].size() -2){ + for ( int i = 0; i < barMap[mapIt.first].size(); i++){ + int coinBar = barMap[mapIt.first][i]; + if (coinBar != iBar && coin_qfineL[i] > 13. && coin_qfineR[i] > 13.){ + float coinEnergy = 0.5 *( energy[coin_chL[i]] + energy[coin_chR[i]]); + if( coinEnergy > energyRanges[label_c][1] && coinEnergy < energyRanges[label_c][2] ){ + std::string label_coin = Form( "1173_bar%02dL-R_vs_1332_bar%02dL-R_%s", iBar, coinBar, mapIt.first.c_str()); + std::string labelTot_coin = Form( "1173_bar%02dL-R_vs_1332_bar%02dL-R_%s_%s", iBar, coinBar, tokens[0].c_str(), tokens[1].c_str()); + std::string label_TProf = Form("bar%02dL-R_%s_%s",iBar, tokens[0].c_str(), tokens[1].c_str()); + + float energyRatioCorr = p1_TProfile_energyRatio_corr_FitFunc[label_TProf][9]->Eval(energyR/energyL) - p1_TProfile_energyRatio_corr_FitFunc[label_TProf][9]->Eval(h1_energyRatio_FitFunc[label_TProf][9]->GetParameter(1)); + h_timeRes_corr_data_coincidence_run[label_coin] -> Fill(timeR -timeL -energyRatioCorr); + h_timeRes_corr_data_coincidence_tot[labelTot_coin] -> Fill(timeR -timeL -energyRatioCorr); + + } + } + } + } + } + + if( avEnergy > energyRanges[label_c][1] && avEnergy < energyRanges[label_c][2] ){ + //coincidence 1332 vs 1173 (no energy in other bars) + int coin_chL[barMap[mapIt.first].size()]; + int coin_chR[barMap[mapIt.first].size()]; + float coin_qfineL[barMap[mapIt.first].size()]; + float coin_qfineR[barMap[mapIt.first].size()]; + + int counter = 0; + for ( int i = 0; i < barMap[mapIt.first].size(); i++){ + int coinBar = barMap[mapIt.first][i]; + coin_chL[i] = chsL[coinBar]; + coin_chR[i] = chsR[coinBar]; + coin_qfineL[i] = qfine[coin_chL[i]]; + coin_qfineR[i] = qfine[coin_chR[i]]; + if ( coin_qfineL[i] < 13. && coin_qfineR[i] < 13.){ + counter ++; + } + } + if(counter == barMap[mapIt.first].size() -2){ + for ( int i = 0; i < barMap[mapIt.first].size(); i++){ + int coinBar = barMap[mapIt.first][i]; + if (coinBar != iBar && coin_qfineL[i] > 13. && coin_qfineR[i] > 13.){ + float coinEnergy = 0.5 *( energy[coin_chL[i]] + energy[coin_chR[i]]); + if( coinEnergy > energyRanges[label_c][0] && coinEnergy < energyRanges[label_c][1] ){ + std::string label_coin = Form( "1332_bar%02dL-R_vs_1173_bar%02dL-R_%s", iBar, coinBar, mapIt.first.c_str()); + std::string labelTot_coin = Form( "1332_bar%02dL-R_vs_1173_bar%02dL-R_%s_%s", iBar, coinBar, tokens[0].c_str(), tokens[1].c_str()); + std::string label_TProf = Form("bar%02dL-R_%s_%s",iBar, tokens[0].c_str(), tokens[1].c_str()); + float energyRatioCorr = p1_TProfile_energyRatio_corr_FitFunc[label_TProf][8]->Eval(energyR/energyL) - p1_TProfile_energyRatio_corr_FitFunc[label_TProf][8]->Eval(h1_energyRatio_FitFunc[label_TProf][8]->GetParameter(1)); + h_timeRes_corr_data_coincidence_run[label_coin] -> Fill(timeR -timeL -energyRatioCorr); + h_timeRes_corr_data_coincidence_tot[labelTot_coin] -> Fill(timeR -timeL -energyRatioCorr); + } + } + } + } + } + } + } + } + + //Draw and save graphs + for ( auto label : c_energy_bar_run){ + label.second -> cd(); + h_energy_bar_run[label.first] -> SetTitle(";energy [a.u.];entries"); + h_energy_bar_run[label.first] -> Draw(); + label.second-> Print(Form("%s/coincidenceResults/energy/runs/c_energy_%s.png",plotDir.c_str(),label.first.c_str())); + label.second-> Print(Form("%s/coincidenceResults/energy/runs/c_energy_%s.pdf",plotDir.c_str(),label.first.c_str())); + outFile -> cd(); + h_energy_bar_run[label.first] -> Write(); + } + for ( auto label : c_energy_bar_tot){ + label.second -> cd(); + h_energy_bar_tot[label.first] -> SetTitle(";energy [a.u.];entries"); + h_energy_bar_tot[label.first] -> Draw(); + label.second-> Print(Form("%s/coincidenceResults/energy/total/c_energy_%s.png",plotDir.c_str(),label.first.c_str())); + label.second-> Print(Form("%s/coincidenceResults/energy/total/c_energy_%s.pdf",plotDir.c_str(),label.first.c_str())); + outFile -> cd(); + h_energy_bar_tot[label.first] -> Write(); + } + for ( auto label : c_energy_1173_bar_run){ + label.second -> cd(); + h_energy_1173_bar_run[label.first] -> SetTitle(";energy [a.u.];entries"); + h_energy_1173_bar_run[label.first] -> Draw(); + label.second-> Print(Form("%s/coincidenceResults/energy/runs/c_1173_peak_%s.png",plotDir.c_str(),label.first.c_str())); + label.second-> Print(Form("%s/coincidenceResults/energy/runs/c_1173_peak_%s.pdf",plotDir.c_str(),label.first.c_str())); + outFile -> cd(); + h_energy_1173_bar_run[label.first] -> Write(); + } + for ( auto label : c_energy_1173_bar_tot){ + label.second -> cd(); + h_energy_1173_bar_tot[label.first] -> SetTitle(";energy [a.u.];entries"); + h_energy_1173_bar_tot[label.first] -> Draw(); + label.second-> Print(Form("%s/coincidenceResults/energy/total/c_1173_peak_%s.png",plotDir.c_str(),label.first.c_str())); + label.second-> Print(Form("%s/coincidenceResults/energy/total/c_1173_peak_%s.pdf",plotDir.c_str(),label.first.c_str())); + outFile -> cd(); + h_energy_1173_bar_tot[label.first] -> Write(); + } + for ( auto label : c_energy_1332_bar_run){ + label.second -> cd(); + h_energy_1332_bar_run[label.first] -> SetTitle(";energy [a.u.];entries"); + h_energy_1332_bar_run[label.first] -> Draw(); + label.second-> Print(Form("%s/coincidenceResults/energy/runs/c_1332_peak_%s.png",plotDir.c_str(),label.first.c_str())); + label.second-> Print(Form("%s/coincidenceResults/energy/runs/c_1332_peak_%s.pdf",plotDir.c_str(),label.first.c_str())); + outFile -> cd(); + h_energy_1332_bar_run[label.first] -> Write(); + } + for ( auto label : c_energy_1332_bar_tot){ + label.second -> cd(); + h_energy_1332_bar_tot[label.first] -> SetTitle(";energy [a.u.];entries"); + h_energy_1332_bar_tot[label.first] -> Draw(); + label.second-> Print(Form("%s/coincidenceResults/energy/total/c_1332_peak_%s.png",plotDir.c_str(),label.first.c_str())); + label.second-> Print(Form("%s/coincidenceResults/energy/total/c_1332_peak_%s.pdf",plotDir.c_str(),label.first.c_str())); + outFile -> cd(); + h_energy_1332_bar_tot[label.first] -> Write(); + } + for ( auto label : c_energy_coincidence_run){ + label.second -> cd(); + h_energy_coincidence_run[label.first] -> SetTitle(";energy [a.u.];entries"); + h_energy_coincidence_run[label.first] -> Draw(); + label.second-> Print(Form("%s/coincidenceResults/energy/runs/c_energy_coincidence_%s.png",plotDir.c_str(),label.first.c_str())); + label.second-> Print(Form("%s/coincidenceResults/energy/runs/c_energy_coincidence_%s.pdf",plotDir.c_str(),label.first.c_str())); + outFile -> cd(); + h_energy_coincidence_run[label.first] -> Write(); + } + for ( auto label : c_energy_coincidence_tot){ + label.second -> cd(); + h_energy_coincidence_tot[label.first] -> SetTitle(";energy [a.u.];entries"); + //h_energy_coincidence_tot[label.first] -> GetYaxis() -> SetRangeUser ( 0, h_energy_coincidence_tot[label.first] -> GetMaximum() *1.1); + h_energy_coincidence_tot[label.first] -> Draw(); + label.second-> Print(Form("%s/coincidenceResults/energy/total/c_energy_coincidence_%s.png",plotDir.c_str(),label.first.c_str())); + label.second-> Print(Form("%s/coincidenceResults/energy/total/c_energy_coincidence_%s.pdf",plotDir.c_str(),label.first.c_str())); + outFile -> cd(); + h_energy_coincidence_tot[label.first] -> Write(); + } + + for ( auto label : c_timeRes_coincidence_run){ + label.second -> cd(); + h_timeRes_coincidence_run[label.first] -> SetTitle(";energy-corrected #Deltat [ps];entries"); + h_timeRes_coincidence_run[label.first] -> Draw(); + label.second-> Print(Form("%s/coincidenceResults/timeResolution/runs/c_energy_coincidence_%s.png",plotDir.c_str(),label.first.c_str())); + label.second-> Print(Form("%s/coincidenceResults/timeResolution/runs/c_energy_coincidence_%s.pdf",plotDir.c_str(),label.first.c_str())); + outFile -> cd(); + h_timeRes_coincidence_run[label.first] -> Write(); + } + + for ( auto label : c_timeRes_coincidence_tot){ + label.second -> cd(); + //raw time resolution + h_timeRes_coincidence_tot[label.first] -> SetTitle(";energy-corrected #Deltat [ps];entries"); + h_timeRes_coincidence_tot[label.first] -> SetLineColor(kRed+1); + float fitXMin = h_timeRes_coincidence_tot[label.first] -> GetBinCenter(h_timeRes_coincidence_tot[label.first]-> GetMaximumBin()) - 200.; + float fitXMax = h_timeRes_coincidence_tot[label.first] -> GetBinCenter(h_timeRes_coincidence_tot[label.first] -> GetMaximumBin()) + 200.; + TF1* fitFunc = new TF1(Form("fitFunc_%s",label.first.c_str()),"gaus",fitXMin,fitXMax); + fitFunc -> SetParameters(1, h_timeRes_coincidence_tot[label.first] -> GetMean(), h_timeRes_coincidence_tot[label.first] -> GetRMS()); + h_timeRes_coincidence_tot[label.first] -> Fit(fitFunc,"QNRSL+",""); + h_timeRes_coincidence_tot[label.first] -> Fit(fitFunc,"QNRSL+","",fitFunc->GetParameter(1)-fitFunc->GetParameter(2),fitFunc->GetParameter(1)+fitFunc->GetParameter(2)); + h_timeRes_coincidence_tot[label.first] -> Draw(); + + float fitXMin2 = fitFunc->GetParameter(1)-fabs(1.5*fitFunc->GetParameter(2)); + float fitXMax2 = fitFunc->GetParameter(1)+fabs(1.5*fitFunc->GetParameter(2)); + TF1* fitFunc2 = new TF1(Form("fitFunc2_energy_%s",label.first.c_str()),"gaus(0)",fitXMin2,fitXMax2); + fitFunc2 -> SetParameter(0,fitFunc->GetParameter(0)); + fitFunc2 -> SetParameter(1,fitFunc->GetParameter(1)); + fitFunc2 -> SetParameter(2,fitFunc->GetParameter(2)); + h_timeRes_coincidence_tot[label.first] -> Fit(fitFunc2,"QRSL+"); + h_timeRes_coincidence_tot[label.first] -> Fit(fitFunc2,"QRSL+"); + + fitFunc2 -> SetLineColor(kRed+1); + fitFunc2 -> SetLineWidth(3); + fitFunc2 -> Draw("same"); + + TLatex* latex = new TLatex(0.20,0.85,Form("{#sigma_{raw}^{gaus} = %.0f ps}",fabs(fitFunc2->GetParameter(2)))); + latex -> SetNDC(); + latex -> SetTextFont(42); + latex -> SetTextSize(0.04); + latex -> SetTextColor(kRed); + latex -> Draw("same"); + + //step 2 corrected time resolution + h_timeRes_corr_coincidence_tot[label.first] -> SetLineColor(kBlue+1); + fitXMin = h_timeRes_corr_coincidence_tot[label.first] -> GetBinCenter(h_timeRes_corr_coincidence_tot[label.first]-> GetMaximumBin()) - 200.; + fitXMax = h_timeRes_corr_coincidence_tot[label.first] -> GetBinCenter(h_timeRes_corr_coincidence_tot[label.first] -> GetMaximumBin()) + 200.; + TF1* fitFunc3 = new TF1(Form("fitFunc_%s",label.first.c_str()),"gaus",fitXMin,fitXMax); + fitFunc3 -> SetParameters(1, h_timeRes_corr_coincidence_tot[label.first] -> GetMean(), h_timeRes_corr_coincidence_tot[label.first] -> GetRMS()); + h_timeRes_corr_coincidence_tot[label.first] -> Fit(fitFunc3,"QNRSL+",""); + h_timeRes_corr_coincidence_tot[label.first] -> Fit(fitFunc3,"QNRSL+","",fitFunc3->GetParameter(1)-fitFunc3->GetParameter(2),fitFunc3->GetParameter(1)+fitFunc3->GetParameter(2)); + h_timeRes_corr_coincidence_tot[label.first] -> Draw("same"); + + fitXMin2 = fitFunc3->GetParameter(1)-fabs(1.5*fitFunc3->GetParameter(2)); + fitXMax2 = fitFunc3->GetParameter(1)+fabs(1.5*fitFunc3->GetParameter(2)); + TF1* fitFunc4 = new TF1(Form("fitFunc2_energyCorr_%s",label.first.c_str()),"gaus(0)",fitXMin2,fitXMax2); + fitFunc4 -> SetParameter(0,fitFunc3->GetParameter(0)); + fitFunc4 -> SetParameter(1,fitFunc3->GetParameter(1)); + fitFunc4 -> SetParameter(2,fitFunc3->GetParameter(2)); + h_timeRes_corr_coincidence_tot[label.first] -> Fit(fitFunc4,"QRSL+"); + h_timeRes_corr_coincidence_tot[label.first] -> Fit(fitFunc4,"QRSL+"); + + fitFunc4 -> SetLineColor(kBlue+1); + fitFunc4 -> SetLineWidth(3); + fitFunc4 -> Draw("same"); + + latex = new TLatex(0.55,0.85,Form("{#sigma_{corr. step2}^{gaus} = %.0f ps}",fabs(fitFunc4->GetParameter(2)))); + latex -> SetNDC(); + latex -> SetTextFont(42); + latex -> SetTextSize(0.04); + latex -> SetTextColor(kBlue); + latex -> Draw("same"); + + + + + //data energy ratio corrected time resolution + h_timeRes_corr_data_coincidence_tot[label.first] -> SetLineColor(kBlue+1); + fitXMin = h_timeRes_corr_data_coincidence_tot[label.first] -> GetBinCenter(h_timeRes_corr_data_coincidence_tot[label.first]-> GetMaximumBin()) - 200.; + fitXMax = h_timeRes_corr_data_coincidence_tot[label.first] -> GetBinCenter(h_timeRes_corr_data_coincidence_tot[label.first] -> GetMaximumBin()) + 200.; + TF1* fitFunc5 = new TF1(Form("fitFunc_%s",label.first.c_str()),"gaus",fitXMin,fitXMax); + fitFunc5 -> SetParameters(1, h_timeRes_corr_data_coincidence_tot[label.first] -> GetMean(), h_timeRes_corr_data_coincidence_tot[label.first] -> GetRMS()); + h_timeRes_corr_data_coincidence_tot[label.first] -> Fit(fitFunc5,"QNRSL+",""); + h_timeRes_corr_data_coincidence_tot[label.first] -> Fit(fitFunc5,"QNRSL+","",fitFunc5->GetParameter(1)-fitFunc5->GetParameter(2),fitFunc5->GetParameter(1)+fitFunc5->GetParameter(2)); + h_timeRes_corr_data_coincidence_tot[label.first] -> Draw("same"); + + fitXMin2 = fitFunc5->GetParameter(1)-fabs(1.5*fitFunc5->GetParameter(2)); + fitXMax2 = fitFunc5->GetParameter(1)+fabs(1.5*fitFunc5->GetParameter(2)); + TF1* fitFunc6 = new TF1(Form("fitFunc2_energyCorr_%s",label.first.c_str()),"gaus(0)",fitXMin2,fitXMax2); + fitFunc6 -> SetParameter(0,fitFunc5->GetParameter(0)); + fitFunc6 -> SetParameter(1,fitFunc5->GetParameter(1)); + fitFunc6 -> SetParameter(2,fitFunc5->GetParameter(2)); + h_timeRes_corr_data_coincidence_tot[label.first] -> Fit(fitFunc6,"QRSL+"); + h_timeRes_corr_data_coincidence_tot[label.first] -> Fit(fitFunc6,"QRSL+"); + + fitFunc6 -> SetLineColor(kGreen+3); + fitFunc6 -> SetLineWidth(3); + fitFunc6 -> Draw("same"); + + latex = new TLatex(0.55,0.70,Form("{#sigma_{corr. data}^{gaus} = %.0f ps}",fabs(fitFunc6->GetParameter(2)))); + latex -> SetNDC(); + latex -> SetTextFont(42); + latex -> SetTextSize(0.04); + latex -> SetTextColor(kGreen+3); + latex -> Draw("same"); + + + label.second-> Print(Form("%s/coincidenceResults/timeResolution/total/c_energy_coincidence_%s.png",plotDir.c_str(),label.first.c_str())); + label.second-> Print(Form("%s/coincidenceResults/timeResolution/total/c_energy_coincidence_%s.pdf",plotDir.c_str(),label.first.c_str())); + outFile -> cd(); + h_timeRes_coincidence_tot[label.first] -> Write(); + h_timeRes_corr_coincidence_tot[label.first] -> Write(); + + + + std::vector tokens = GetTokens(label.first,'_'); + int iBin = 0; + std::string label_graph = Form("%s_%s_%s", tokens[1].c_str(), tokens[5].c_str(), tokens[6].c_str()); + if( atof(tokens[0].c_str()) == 1173){ + iBin = 8; + } + if( atof(tokens[0].c_str()) == 1332){ + iBin = 9; + } + TGraphErrors* graph = new TGraphErrors(); + graph -> SetPoint( 0, gr_timeRes_summary[label_graph][iBin].size() +1, fitFunc4->GetParameter(2)); + graph -> SetPointError( 0, 0, 0); + graph -> SetMarkerColor (2* (gr_timeRes_summary[label_graph][iBin].size() +1)); + graph -> SetMarkerStyle(26); + gr_timeRes_summary[label_graph][iBin].push_back( (TGraphErrors*) graph); + l_timeRes_summary[label_graph][iBin] -> AddEntry( graph, Form("%s_coincidence_step2Corr",tokens[4].c_str()), "p"); + + TGraphErrors* graph2 = new TGraphErrors(); + graph2 -> SetPoint( 0, gr_timeRes_summary[label_graph][iBin].size() +1, fitFunc6->GetParameter(2)); + graph2 -> SetPointError( 0, 0, 0); + graph2 -> SetMarkerColor (2*(gr_timeRes_summary[label_graph][iBin].size() +1)-1); + graph2 -> SetMarkerStyle(24); + gr_timeRes_summary[label_graph][iBin].push_back( (TGraphErrors*) graph2); + l_timeRes_summary[label_graph][iBin] -> AddEntry( graph2, Form("%s_coincidence_dataCorr",tokens[4].c_str()), "p"); + + } + + for( auto label: c_step2_TProfile_run){ + label.second -> cd(); + label.second -> Print(Form("%s/coincidenceResults/timeResolution/ratioEnergyCorrection/c_deltaT_vs_energyRatio_%s.png",plotDir.c_str(),label.first.c_str())); + label.second -> Print(Form("%s/coincidenceResults/timeResolution/ratioEnergyCorrection/c_deltaT_vs_energyRatio_%s.pdf",plotDir.c_str(),label.first.c_str())); + } + + for( auto label: c_TProfile_energyRatio_corr){ + for( int iBin = 8; iBin < 10; iBin++){ + + label.second[iBin] -> cd(); + label.second[iBin] -> Print(Form("%s/coincidenceResults/timeResolution/ratioEnergyCorrection/c_deltaT_vs_energyRatio_%s_energyBin%d_coincidenceData.png",plotDir.c_str(),label.first.c_str(),iBin)); + label.second[iBin] -> Print(Form("%s/coincidenceResults/timeResolution/ratioEnergyCorrection/c_deltaT_vs_energyRatio_%s_energyBin%d_coincidenceData.pdf",plotDir.c_str(),label.first.c_str(),iBin)); + } + } + + + for ( auto label: c_timeRes_summary){ + for( int iBin = 8; iBin < 10; iBin++){ + label.second[iBin] -> cd(); + gr_timeRes_summary[label.first][iBin][gr_timeRes_summary[label.first][iBin].size()-1] -> SetTitle(";;energy-corrected #Deltat [ps]"); + gr_timeRes_summary[label.first][iBin][gr_timeRes_summary[label.first][iBin].size()-1] -> GetHistogram() -> GetXaxis() -> SetLimits(0,gr_timeRes_summary[label.first][iBin].size() + 10); + gr_timeRes_summary[label.first][iBin][gr_timeRes_summary[label.first][iBin].size()-1] -> GetYaxis() -> SetRangeUser(0, 600); + gr_timeRes_summary[label.first][iBin][gr_timeRes_summary[label.first][iBin].size()-1] -> Draw("AP"); + for ( int i=0; i < gr_timeRes_summary[label.first][iBin].size()-1; i++){ + gr_timeRes_summary[label.first][iBin][i] -> Draw("psame"); + } + l_timeRes_summary[label.first][iBin]-> Draw("same"); + label.second[iBin] -> Print(Form("%s/coincidenceResults/timeResolution/timeResolutionComparison/c_timeResolutioComparison_%s_energyBin%d.png",plotDir.c_str(),label.first.c_str(),iBin)); + label.second[iBin] -> Print(Form("%s/coincidenceResults/timeResolution/timeResolutionComparison/c_timeResolutioComparison_%s_energyBin%d.pdf",plotDir.c_str(),label.first.c_str(),iBin)); + } + } +} + + diff --git a/main/moduleCharacterization_step1.cpp b/main/moduleCharacterization_step1.cpp index 293fb479..9b162089 100644 --- a/main/moduleCharacterization_step1.cpp +++ b/main/moduleCharacterization_step1.cpp @@ -45,7 +45,6 @@ int main(int argc, char** argv) //--- parse the config file CfgManager opts; opts.ParseConfigFile(argv[1]); - int debugMode = 0; if( argc > 2 ) debugMode = atoi(argv[2]); @@ -58,6 +57,7 @@ int main(int argc, char** argv) int maxEntries = opts.GetOpt("Input.maxEntries"); int usePedestals = opts.GetOpt("Input.usePedestals"); TChain* tree = new TChain("data","data"); + std::stringstream ss(runs); std::string token; @@ -147,13 +147,25 @@ int main(int argc, char** argv) //--- get plot settings std::vector Vov = opts.GetOpt >("Plots.Vov"); - std::vector energyBins = opts.GetOpt >("Plots.energyBins"); + std::vector energyBins; + if(!opts.GetOpt("Input.sourceName").compare("Na22")){ + energyBins = opts.GetOpt >("Plots.energyBinsNa22"); + } + if(!opts.GetOpt("Input.sourceName").compare("Co60")){ + energyBins = opts.GetOpt >("Plots.energyBinsCo60"); + } std::map map_energyBins; for(unsigned int ii = 0; ii < Vov.size(); ++ii) map_energyBins[Vov[ii]] = energyBins[ii]; float energyMin = opts.GetOpt("Plots.energyMin"); - float energyMax = opts.GetOpt("Plots.energyMax"); + float energyMax; + if(!opts.GetOpt("Input.sourceName").compare("Na22")){ + energyMax = opts.GetOpt("Plots.energyMaxNa22"); + } + if(!opts.GetOpt("Input.sourceName").compare("Co60")){ + energyMax = opts.GetOpt("Plots.energyMaxCo60"); + } @@ -176,95 +188,229 @@ int main(int argc, char** argv) //------------------------ //--- 1st loop over events - ModuleEventClass anEvent; - - int nEntries = tree->GetEntries(); - if( maxEntries > 0 ) nEntries = maxEntries; - for(int entry = 0; entry < nEntries; ++entry) - { - tree -> GetEntry(entry); - if( entry%200000 == 0 ) - { - std::cout << "\n>>> 1st loop: reading entry " << entry << " / " << nEntries << " (" << 100.*entry/nEntries << "%)" << std::endl; - TrackProcess(cpu, mem, vsz, rss); - } - - float Vov = step1; - float vth1 = float(int(step2/10000)-1);; - // float vth2 = float(int((step2-10000*vth1)/100)-1); - // float vthe = float(int((step2-10000*vth1-step2-100*vth2)/1)-1); - - - for(int iBar = 0; iBar < 16; ++iBar) - { - int index( (10000*int(Vov*100.)) + (100*vth1) + iBar ); - - - int chL(chsL[iBar]); - int chR(chsR[iBar]); - - float qfineL(qfine[chL]); - float qfineR(qfine[chR]); - float totL(0.001*tot[chL]); - float totR(0.001*tot[chR]); - - if( totL <= 0. || totR <= 0. ) continue; - if( qfineL < 13. || qfineR < 13. ) continue; - - float energyL(energy[chL]); - float energyR(energy[chR]); - long long timeL(time[chL]); - long long timeR(time[chR]); - - - //--- create histograms, if needed - if( h1_totL[index] == NULL ) - { - h1_qfineL[index] = new TH1F(Form("h1_qfine_bar%02dL_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",512,-0.5,511.5); - h1_qfineR[index] = new TH1F(Form("h1_qfine_bar%02dR_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",512,-0.5,511.5); - - h1_totL[index] = new TH1F(Form("h1_tot_bar%02dL_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",2000,0.,1000.); - h1_totR[index] = new TH1F(Form("h1_tot_bar%02dR_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",2000,0.,1000.); - - h1_energyL[index] = new TH1F(Form("h1_energy_bar%02dL_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",map_energyBins[Vov],energyMin,energyMax); - h1_energyR[index] = new TH1F(Form("h1_energy_bar%02dR_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",map_energyBins[Vov],energyMin,energyMax); - - outTrees[index] = new TTree(Form("data_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1),Form("data_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1)); - outTrees[index] -> Branch("event",&anEvent); - - h1_energyLR[index] = new TH1F(Form("h1_energy_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",map_energyBins[Vov],energyMin,energyMax); - } - - - //--- fill histograms - h1_qfineL[index] -> Fill( qfineL ); - h1_totL[index] -> Fill( totL ); - h1_energyL[index] -> Fill( energyL ); - - h1_qfineR[index] -> Fill( qfineR ); - h1_totR[index] -> Fill( totR ); - h1_energyR[index] -> Fill( energyR ); - - h1_energyLR[index] -> Fill(0.5*(energyL+energyR)); - - anEvent.barID = iBar; - anEvent.Vov = Vov; - anEvent.vth1 = vth1; - anEvent.energyL = energyL; - anEvent.energyR = energyR; - anEvent.timeL = timeL; - anEvent.timeR = timeR; - outTrees[index] -> Fill(); - } - } - std::cout << std::endl; - - - - int bytes = outFile -> Write(); - std::cout << "============================================" << std::endl; - std::cout << "nr of B written: " << int(bytes) << std::endl; - std::cout << "nr of KB written: " << int(bytes/1024.) << std::endl; - std::cout << "nr of MB written: " << int(bytes/1024./1024.) << std::endl; - std::cout << "============================================" << std::endl; + + if (!opts.GetOpt("Input.sourceName").compare("Na22")){ + + ModuleEventClass anEvent; + + float qfineL[16]; + float qfineR[16]; + float totL[16]; + float totR[16]; + long long int timeL[16]; + long long int timeR[16]; + float energyL[16]; + float energyR[16]; + int chL[16]; + int chR[16]; + + int nEntries = tree->GetEntries(); + if( maxEntries > 0 ) nEntries = maxEntries; + for(int entry = 0; entry < nEntries; ++entry) + { + tree -> GetEntry(entry); + if( entry%200000 == 0 ) + { + std::cout << "\n>>> 1st loop: reading entry " << entry << " / " << nEntries << " (" << 100.*entry/nEntries << "%)" << std::endl; + TrackProcess(cpu, mem, vsz, rss); + } + + float Vov = step1; + float vth1 = float(int(step2/10000)-1);; + // float vth2 = float(int((step2-10000*vth1)/100)-1); + // float vthe = float(int((step2-10000*vth1-step2-100*vth2)/1)-1); + + + for(int iBar = 0; iBar < 16; ++iBar) + { + + + + chL[iBar]=chsL[iBar]; + chR[iBar]=chsR[iBar]; + + + + + qfineL[iBar]=qfine[iBar]; + qfineR[iBar]=qfine[iBar]; + totL[iBar]=0.001*tot[iBar]; + totR[iBar]=0.001*tot[iBar]; + + energyL[iBar]=energy[iBar]; + energyR[iBar]=energy[iBar]; + timeL[iBar]=time[iBar]; + timeR[iBar]=time[iBar]; + + } + + int maxEn=0; + int maxBar=0; + + + for(int iBar = 0; iBar < 16; ++iBar) + { + + if (qfineL[iBar]>13 && qfineR[iBar]>13){ + float energyMean=(energyL[iBar]+energyR[iBar])/2; + if(energyMean>maxEn){ + maxEn = energyMean; + maxBar = iBar; + } + } + + int index( (10000*int(Vov*100.)) + (100*vth1) + iBar ); + + //--- create histograms, if needed + if( h1_totL[index] == NULL ) + { + h1_qfineL[index] = new TH1F(Form("h1_qfine_bar%02dL_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",512,-0.5,511.5); + h1_qfineR[index] = new TH1F(Form("h1_qfine_bar%02dR_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",512,-0.5,511.5); + + h1_totL[index] = new TH1F(Form("h1_tot_bar%02dL_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",2000,0.,1000.); + h1_totR[index] = new TH1F(Form("h1_tot_bar%02dR_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",2000,0.,1000.); + + h1_energyL[index] = new TH1F(Form("h1_energy_bar%02dL_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",map_energyBins[Vov],energyMin,energyMax); + h1_energyR[index] = new TH1F(Form("h1_energy_bar%02dR_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",map_energyBins[Vov],energyMin,energyMax); + + outTrees[index] = new TTree(Form("data_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1),Form("data_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1)); + outTrees[index] -> Branch("event",&anEvent); + + h1_energyLR[index] = new TH1F(Form("h1_energy_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",map_energyBins[Vov],energyMin,energyMax); + } + + } + int index( (10000*int(Vov*100.)) + (100*vth1) + maxBar ); + + //--- fill histograms + h1_qfineL[index] -> Fill( qfineL[maxBar] ); + h1_totL[index] -> Fill( totL[maxBar] ); + h1_energyL[index] -> Fill( energyL[maxBar] ); + + h1_qfineR[index] -> Fill( qfineR[maxBar] ); + h1_totR[index] -> Fill( totR[maxBar] ); + h1_energyR[index] -> Fill( energyR[maxBar] ); + + h1_energyLR[index] -> Fill(0.5*(energyL[maxBar]+energyR[maxBar])); + + anEvent.barID = maxBar; + anEvent.Vov = Vov; + anEvent.vth1 = vth1; + anEvent.energyL = energyL[maxBar]; + anEvent.energyR = energyR[maxBar]; + anEvent.timeL = timeL[maxBar]; + anEvent.timeR = timeR[maxBar]; + outTrees[index] -> Fill(); + + } + std::cout << std::endl; + + + + int bytes = outFile -> Write(); + std::cout << "============================================" << std::endl; + std::cout << "nr of B written: " << int(bytes) << std::endl; + std::cout << "nr of KB written: " << int(bytes/1024.) << std::endl; + std::cout << "nr of MB written: " << int(bytes/1024./1024.) << std::endl; + std::cout << "============================================" << std::endl; + + } + + if (!opts.GetOpt("Input.sourceName").compare("Co60")){ + ModuleEventClass anEvent; + + int nEntries = tree->GetEntries(); + if( maxEntries > 0 ) nEntries = maxEntries; + for(int entry = 0; entry < nEntries; ++entry) + { + tree -> GetEntry(entry); + if( entry%200000 == 0 ) + { + std::cout << "\n>>> 1st loop: reading entry " << entry << " / " << nEntries << " (" << 100.*entry/nEntries << "%)" << std::endl; + TrackProcess(cpu, mem, vsz, rss); + } + + float Vov = step1; + float vth1 = float(int(step2/10000)-1);; + // float vth2 = float(int((step2-10000*vth1)/100)-1); + // float vthe = float(int((step2-10000*vth1-step2-100*vth2)/1)-1); + + + for(int iBar = 0; iBar < 16; ++iBar) + { + int index( (10000*int(Vov*100.)) + (100*vth1) + iBar ); + + + int chL(chsL[iBar]); + int chR(chsR[iBar]); + + float qfineL(qfine[chL]); + float qfineR(qfine[chR]); + float totL(0.001*tot[chL]); + float totR(0.001*tot[chR]); + + if( totL <= 0. || totR <= 0. ) continue; + if( qfineL < 13. || qfineR < 13. ) continue; + + float energyL(energy[chL]); + float energyR(energy[chR]); + long long timeL(time[chL]); + long long timeR(time[chR]); + + + //--- create histograms, if needed + if( h1_totL[index] == NULL ) + { + h1_qfineL[index] = new TH1F(Form("h1_qfine_bar%02dL_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",512,-0.5,511.5); + h1_qfineR[index] = new TH1F(Form("h1_qfine_bar%02dR_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",512,-0.5,511.5); + + h1_totL[index] = new TH1F(Form("h1_tot_bar%02dL_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",2000,0.,1000.); + h1_totR[index] = new TH1F(Form("h1_tot_bar%02dR_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",2000,0.,1000.); + + h1_energyL[index] = new TH1F(Form("h1_energy_bar%02dL_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",map_energyBins[Vov],energyMin,energyMax); + h1_energyR[index] = new TH1F(Form("h1_energy_bar%02dR_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",map_energyBins[Vov],energyMin,energyMax); + + outTrees[index] = new TTree(Form("data_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1),Form("data_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1)); + outTrees[index] -> Branch("event",&anEvent); + + h1_energyLR[index] = new TH1F(Form("h1_energy_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",map_energyBins[Vov],energyMin,energyMax); + } + + + //--- fill histograms + h1_qfineL[index] -> Fill( qfineL ); + h1_totL[index] -> Fill( totL ); + h1_energyL[index] -> Fill( energyL ); + + h1_qfineR[index] -> Fill( qfineR ); + h1_totR[index] -> Fill( totR ); + h1_energyR[index] -> Fill( energyR ); + + h1_energyLR[index] -> Fill(0.5*(energyL+energyR)); + + anEvent.barID = iBar; + anEvent.Vov = Vov; + anEvent.vth1 = vth1; + anEvent.energyL = energyL; + anEvent.energyR = energyR; + anEvent.timeL = timeL; + anEvent.timeR = timeR; + outTrees[index] -> Fill(); + } + } + std::cout << std::endl; + + + + int bytes = outFile -> Write(); + std::cout << "============================================" << std::endl; + std::cout << "nr of B written: " << int(bytes) << std::endl; + std::cout << "nr of KB written: " << int(bytes/1024.) << std::endl; + std::cout << "nr of MB written: " << int(bytes/1024./1024.) << std::endl; + std::cout << "============================================" << std::endl; + } + + + + } diff --git a/main/moduleCharacterization_step2.cpp b/main/moduleCharacterization_step2.cpp index 33c12e47..42315710 100644 --- a/main/moduleCharacterization_step2.cpp +++ b/main/moduleCharacterization_step2.cpp @@ -158,6 +158,7 @@ int main(int argc, char** argv) stepLabels.erase(std::unique(stepLabels.begin(),stepLabels.end()),stepLabels.end()); + //--- define histograms std::string outFileName = opts.GetOpt("Output.outFileNameStep2"); TFile* outFile = TFile::Open(outFileName.c_str(),"RECREATE"); @@ -173,6 +174,7 @@ int main(int argc, char** argv) float energy1275L; float theIndex; float timeRes; + float energyBinValCo60[3]; std::map h1_energyRatio; @@ -215,8 +217,7 @@ int main(int argc, char** argv) std::string source = opts.GetOpt("Input.sourceName"); std::string Na22 = "Na22"; std::string Co60 = "Co60"; - for(auto stepLabel : stepLabels) - { + for(auto stepLabel : stepLabels){ float Vov = map_Vovs[stepLabel]; float vth1 = map_ths[stepLabel]; std::string VovLabel(Form("Vov%.1f",Vov)); @@ -226,9 +227,9 @@ int main(int argc, char** argv) //-------------------------------------------------------- //loop sulle barre - for(int iBar = 0; iBar < 16; ++iBar) - { + for(int iBar = 0; iBar < 16; ++iBar){ int index( (10000*int(Vov*100.)) + (100*vth1) + iBar ); + outTrees[index] = new TTree(Form("data_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1),Form("data_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1)); outTrees[index] -> Branch("energyPeak511LR",&energy511LR); @@ -238,22 +239,23 @@ int main(int argc, char** argv) outTrees[index] -> Branch("energyPeak511R",&energy511R); outTrees[index] -> Branch("energyPeak1275R",&energy1275R); outTrees[index] -> Branch("indexID",&theIndex); - - - for(auto LRLabel : LRLabels ) - { + + if(!source.compare(Co60)){ + outTrees[index] -> Branch("energyBinValCo60",energyBinValCo60, "energyBinValCo60[3]/F"); + } + for(auto LRLabel : LRLabels ){ //label histo - std::string label(Form("bar%02d%s_%s",iBar,LRLabel.c_str(),stepLabel.c_str())); + std::string label(Form("bar%02d%s_%s",iBar,LRLabel.c_str(),stepLabel.c_str())); histo = (TH1F*)( inFile->Get(Form("h1_qfine_%s",label.c_str())) ); if( !histo ) continue; if( histo->GetEntries() < 100 ) continue; - latex = new TLatex(0.40,0.85,Form("#splitline{bar %02d%s}{V_{OV} = %.1f V, th. = %d DAC}",iBar,LRLabel.c_str(),Vov,int(vth1))); + latex = new TLatex(0.50,0.90,Form("#splitline{bar %02d%s}{V_{OV} = %.1f V, th. = %d DAC}",iBar,LRLabel.c_str(),Vov,int(vth1))); latex -> SetNDC(); latex -> SetTextFont(42); - latex -> SetTextSize(0.04); + latex -> SetTextSize(0.035); latex -> SetTextColor(kRed); @@ -301,213 +303,217 @@ int main(int argc, char** argv) histo -> SetLineWidth(2); histo -> Draw(); - - if( LRLabel.compare("L")){ - if(!source.compare(Na22)){ - rangesL[index] = new std::vector; - peaksL[index] = Na22SpectrumAnalyzer(histo,rangesL[index]); - - if (peaksL[index]["0.511 MeV"].first== -9999){ - peaksL.erase(index); - rangesL.erase(index); - peaksL[index]["0.511 MeV"].first = -10; - peaksL[index]["1.275 MeV"].first = -10; - } + if( LRLabel.compare("L")){ + if(!source.compare(Na22)){ + rangesL[index] = new std::vector; + peaksL[index] = Na22SpectrumAnalyzer(histo,rangesL[index]); + + if (peaksL[index]["0.511 MeV"].first > 0){ + for(auto peak : peaksL[index] ) + { + histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first); + break; + } + } + if (peaksL[index]["0.511 MeV"].first== -9999){ + for(auto peak : peaksL[index] ) + { + histo -> GetXaxis() -> SetRangeUser(0.,40); + break; + } + peaksL.erase(index); + rangesL.erase(index); + peaksL[index]["0.511 MeV"].first = -10; + peaksL[index]["1.275 MeV"].first = -10; + + } - energy511L = peaksL[index]["0.511 MeV"].first; - energy1275L = peaksL[index]["1.275 MeV"].first; - theIndex = index; - outTrees[index] -> Fill(); - - - - for(auto peak : peaksL[index] ) - { - histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first); - break; - } - - - histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum()); - if(peaksL[index]["0.511 MeV"].first != -10){ - for(int i = 1; i < rangesL[index]->size(); i++){ - TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesL[index]->at(i)) - histo->FindBin(rangesL[index]-> at(i-1)), rangesL[index]-> at(i-1), rangesL[index]->at(i)); - int j = 1; - for (int bin = histo->FindBin(rangesL[index]->at(i-1)) ; bin < histo -> FindBin(rangesL[index]->at(i)) +1 ; bin++){ - binHisto -> SetBinContent( j, histo->GetBinContent(bin)); - j++; + energy511L = peaksL[index]["0.511 MeV"].first; + energy1275L = peaksL[index]["1.275 MeV"].first; + theIndex = index; + outTrees[index] -> Fill(); + + + histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum()); + if(peaksL[index]["0.511 MeV"].first != -10){ + for(int i = 1; i < rangesL[index]->size(); i++){ + TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesL[index]->at(i)) - histo->FindBin(rangesL[index]-> at(i-1)), rangesL[index]-> at(i-1), rangesL[index]->at(i)); + int j = 1; + for (int bin = histo->FindBin(rangesL[index]->at(i-1)) ; bin < histo -> FindBin(rangesL[index]->at(i)) +1 ; bin++){ + binHisto -> SetBinContent( j, histo->GetBinContent(bin)); + j++; + } + energyBinL[index][i] = binHisto -> GetMean(); + binHisto -> Delete(); + } + } } - energyBinL[index][i] = binHisto -> GetMean(); - binHisto -> Delete(); - } - } - - } + if(!source.compare(Co60)){ + rangesL[index] = new std::vector; + peaksL[index] = Co60SpectrumAnalyzer(histo,rangesL[index]); + + if (peaksL[index]["1.173 MeV"].first > 0){ + for(auto peak : peaksL[index] ) + { + histo -> GetXaxis() -> SetRangeUser(0.,2.*peak.second.first); + break; + } + } - if(!source.compare(Co60)){ - rangesL[index] = new std::vector; - peaksL[index] = Co60SpectrumAnalyzer(histo,rangesL[index]); - - - if (peaksL[index]["1.173 MeV"].first== -9999){ - peaksL.erase(index); - rangesL.erase(index); - peaksL[index]["1.173 MeV"].first = -10; - peaksL[index]["1.332 MeV"].first = -10; - - } + if (peaksL[index]["1.173 MeV"].first== -9999){ + for(auto peak : peaksL[index] ) + { + histo -> GetXaxis() -> SetRangeUser(0.,50.); + break; + } + peaksL.erase(index); + rangesL.erase(index); + peaksL[index]["1.173 MeV"].first = -10; + peaksL[index]["1.332 MeV"].first = -10; + + } + energy511L = peaksL[index]["1.173 MeV"].first; + energy1275L = peaksL[index]["1.332 MeV"].first; + theIndex = index; + outTrees[index] -> Fill(); + - energy511L = peaksL[index]["1.173 MeV"].first; - energy1275L = peaksL[index]["1.332 MeV"].first; - theIndex = index; - outTrees[index] -> Fill(); - - - - for(auto peak : peaksL[index] ) - { - histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first); - break; - } - - - histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum()); - if(peaksL[index]["1.173 MeV"].first != -10){ - for(int i = 1; i < rangesL[index]->size(); i++){ - TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesL[index]->at(i)) - histo->FindBin(rangesL[index]-> at(i-1)), rangesL[index]-> at(i-1), rangesL[index]->at(i)); - int j = 1; - for (int bin = histo->FindBin(rangesL[index]->at(i-1)) ; bin < histo -> FindBin(rangesL[index]->at(i)) +1 ; bin++){ - binHisto -> SetBinContent( j, histo->GetBinContent(bin)); - j++; + + if(peaksL[index]["1.173 MeV"].first != -10){ + for(int i = 1; i < rangesL[index]->size(); i++){ + TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesL[index]->at(i)) - histo->FindBin(rangesL[index]-> at(i-1)), rangesL[index]-> at(i-1), rangesL[index]->at(i)); + int j = 1; + for (int bin = histo->FindBin(rangesL[index]->at(i-1)) ; bin < histo -> FindBin(rangesL[index]->at(i)) +1 ; bin++){ + binHisto -> SetBinContent( j, histo->GetBinContent(bin)); + j++; + } + energyBinL[index][i] = binHisto -> GetMean(); + binHisto -> Delete(); + } + } + + } + if( source.compare(Na22) && source.compare(Co60)){ + std::cout << " Missspelled radioactive source " << std::endl; + return(0); } - energyBinL[index][i] = binHisto -> GetMean(); - binHisto -> Delete(); } - } - - } - - - if( source.compare(Na22) && source.compare(Co60)){ - std::cout << " Missspelled radioactive source " << std::endl; - return(0); - } - } + + if( LRLabel.compare("R")){ + if(!source.compare(Na22)){ + + rangesR[index] = new std::vector; + peaksR[index] = Na22SpectrumAnalyzer(histo,rangesR[index]); + + if (peaksR[index]["0.511 MeV"].first > 0){ + for(auto peak : peaksR[index] ) + { + histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first); + break; + } + } - if( LRLabel.compare("R")){ - if(!source.compare(Na22)){ - rangesR[index] = new std::vector; - peaksR[index] = Na22SpectrumAnalyzer(histo,rangesR[index]); + if (peaksR[index]["0.511 MeV"].first== -9999){ + for(auto peak : peaksR[index] ) + { + histo -> GetXaxis() -> SetRangeUser(0.,40.); + break; + } + peaksR.erase(index); + rangesR.erase(index); + peaksR[index]["0.511 MeV"].first = -10; + peaksR[index]["1.275 MeV"].first = -10; + } - if (peaksR[index]["0.511 MeV"].first== -9999){ - peaksR.erase(index); - rangesR.erase(index); - peaksR[index]["0.511 MeV"].first = -10; - peaksR[index]["1.275 MeV"].first = -10; - - } + energy511R = peaksR[index]["0.511 MeV"].first; + energy1275R = peaksR[index]["1.275 MeV"].first; + theIndex = index; + outTrees[index] -> Fill(); + + + + histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum()); + if(peaksR[index]["0.511 MeV"].first != -10){ + for(int i = 1; i < rangesR[index]->size(); i++){ + TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesR[index]->at(i)) - histo->FindBin(rangesR[index]-> at(i-1)), rangesR[index]-> at(i-1), rangesR[index]->at(i)); + int j = 1; + for (int bin = histo->FindBin(rangesR[index]->at(i-1)) ; bin < histo -> FindBin(rangesR[index]->at(i)) +1 ; bin++){ + binHisto -> SetBinContent( j, histo->GetBinContent(bin)); + j++; + } + energyBinR[index][i] = binHisto -> GetMean(); + binHisto -> Delete(); + } + } + } + + if(!source.compare(Co60)){ + rangesR[index] = new std::vector; + peaksR[index] = Co60SpectrumAnalyzer(histo,rangesR[index]); + + if (peaksR[index]["1.173 MeV"].first > 0){ + for(auto peak : peaksR[index] ) + { + histo -> GetXaxis() -> SetRangeUser(0.,2.*peak.second.first); + break; + } + } + if (peaksR[index]["1.173 MeV"].first== -9999){ + for(auto peak : peaksR[index] ) + { + histo -> GetXaxis() -> SetRangeUser(0.,50.); + break; + } + peaksR.erase(index); + rangesR.erase(index); + peaksR[index]["1.173 MeV"].first = -10; + peaksR[index]["1.332 MeV"].first = -10; + + } - energy511R = peaksR[index]["0.511 MeV"].first; - energy1275R = peaksR[index]["1.275 MeV"].first; - theIndex = index; - outTrees[index] -> Fill(); - - - - for(auto peak : peaksR[index] ) - { - histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first); - break; - } - - - histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum()); - if(peaksR[index]["0.511 MeV"].first != -10){ - for(int i = 1; i < rangesR[index]->size(); i++){ - TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesR[index]->at(i)) - histo->FindBin(rangesR[index]-> at(i-1)), rangesR[index]-> at(i-1), rangesR[index]->at(i)); - int j = 1; - for (int bin = histo->FindBin(rangesR[index]->at(i-1)) ; bin < histo -> FindBin(rangesR[index]->at(i)) +1 ; bin++){ - binHisto -> SetBinContent( j, histo->GetBinContent(bin)); - j++; + + energy511R = peaksR[index]["1.173 MeV"].first; + energy1275R = peaksR[index]["1.332 MeV"].first; + theIndex = index; + outTrees[index] -> Fill(); + + + + if(peaksR[index]["1.173 MeV"].first != -10){ + for(int i = 1; i < rangesR[index]->size(); i++){ + TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesR[index]->at(i)) - histo->FindBin(rangesR[index]-> at(i-1)), rangesR[index]-> at(i-1), rangesR[index]->at(i)); + int j = 1; + for (int bin = histo->FindBin(rangesR[index]->at(i-1)) ; bin < histo -> FindBin(rangesR[index]->at(i)) +1 ; bin++){ + binHisto -> SetBinContent( j, histo->GetBinContent(bin)); + j++; + } + energyBinR[index][i] = binHisto -> GetMean(); + binHisto -> Delete(); + } + } } - energyBinR[index][i] = binHisto -> GetMean(); - binHisto -> Delete(); } - } + latex -> Draw("same"); + histo -> Write(); + c -> Print(Form("%s/energy/c_energy__%s.png",plotDir.c_str(),label.c_str())); + c -> Print(Form("%s/energy/c_energy__%s.pdf",plotDir.c_str(),label.c_str())); + delete c; } - if(!source.compare(Co60)){ - rangesR[index] = new std::vector; - peaksR[index] = Co60SpectrumAnalyzer(histo,rangesR[index]); - - - if (peaksR[index]["1.173 MeV"].first== -9999){ - peaksR.erase(index); - rangesR.erase(index); - peaksR[index]["1.173 MeV"].first = -10; - peaksR[index]["1.332 MeV"].first = -10; - - } - - - energy511R = peaksR[index]["1.173 MeV"].first; - energy1275R = peaksR[index]["1.332 MeV"].first; - theIndex = index; - outTrees[index] -> Fill(); - - - - for(auto peak : peaksR[index] ) + { - histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first); - break; - } - - - histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum()); - if(peaksR[index]["1.173 MeV"].first != -10){ - for(int i = 1; i < rangesR[index]->size(); i++){ - TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesR[index]->at(i)) - histo->FindBin(rangesR[index]-> at(i-1)), rangesR[index]-> at(i-1), rangesR[index]->at(i)); - int j = 1; - for (int bin = histo->FindBin(rangesR[index]->at(i-1)) ; bin < histo -> FindBin(rangesR[index]->at(i)) +1 ; bin++){ - binHisto -> SetBinContent( j, histo->GetBinContent(bin)); - j++; - } - energyBinR[index][i] = binHisto -> GetMean(); - binHisto -> Delete(); - } - } - - } - - - if( source.compare(Na22) && source.compare(Co60)){ - std::cout << " Missspelled radioactive source " << std::endl; - return(0); - } - } - - - latex -> Draw("same"); - histo -> Write(); - c -> Print(Form("%s/energy/c_energy__%s.png",plotDir.c_str(),label.c_str())); - c -> Print(Form("%s/energy/c_energy__%s.pdf",plotDir.c_str(),label.c_str())); - delete c; - } - - - { std::string label(Form("bar%02dL-R_%s",iBar,stepLabel.c_str())); + - latex = new TLatex(0.40,0.85,Form("#splitline{bar %02d}{V_{OV} = %.1f V, th. = %d DAC}",iBar,Vov,int(vth1))); + latex = new TLatex(0.50,0.90,Form("#splitline{bar %02d}{V_{OV} = %.1f V, th. = %d DAC}",iBar,Vov,int(vth1))); latex -> SetNDC(); latex -> SetTextFont(42); - latex -> SetTextSize(0.04); + latex -> SetTextSize(0.035); latex -> SetTextColor(kRed); @@ -522,116 +528,120 @@ int main(int argc, char** argv) histo -> SetLineWidth(2); histo -> SetLineColor(kRed); histo -> Draw(); - - + + - if(!source.compare(Na22)){ - rangesLR[index] = new std::vector; - peaksLR[index] = Na22SpectrumAnalyzer(histo,rangesLR[index]); + if(!source.compare(Na22)){ + rangesLR[index] = new std::vector; + peaksLR[index] = Na22SpectrumAnalyzer(histo,rangesLR[index]); + + if (peaksLR[index]["0.511 MeV"].first > 0){ + for(auto peak : peaksLR[index] ) + { + histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first); + break; + } + } - if (peaksLR[index]["0.511 MeV"].first== -9999){ - peaksLR.erase(index); - rangesLR.erase(index); - peaksLR[index]["0.511 MeV"].first = -10; - peaksLR[index]["1.275 MeV"].first = -10; - - } + if (peaksLR[index]["0.511 MeV"].first== -9999){ + for(auto peak : peaksLR[index] ) + { + histo -> GetXaxis() -> SetRangeUser(0.,40.); + break; + } + peaksLR.erase(index); + rangesLR.erase(index); + peaksLR[index]["0.511 MeV"].first = -10; + peaksLR[index]["1.275 MeV"].first = -10; + + } - energy511LR = peaksLR[index]["0.511 MeV"].first; - energy1275LR = peaksLR[index]["1.275 MeV"].first; - theIndex = index; - outTrees[index] -> Fill(); - - - - for(auto peak : peaksLR[index] ) - { - histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first); - break; - } - - - histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum()); - if(peaksLR[index]["0.511 MeV"].first != -10){ - for(int i = 1; i < rangesLR[index]->size(); i++){ - TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesLR[index]->at(i)) - histo->FindBin(rangesLR[index]-> at(i-1)), rangesLR[index]-> at(i-1), rangesLR[index]->at(i)); - int j = 1; - for (int bin = histo->FindBin(rangesLR[index]->at(i-1)) ; bin < histo -> FindBin(rangesLR[index]->at(i)) +1 ; bin++){ - binHisto -> SetBinContent( j, histo->GetBinContent(bin)); - j++; + energy511LR = peaksLR[index]["0.511 MeV"].first; + energy1275LR = peaksLR[index]["1.275 MeV"].first; + theIndex = index; + outTrees[index] -> Fill(); + + + + + if(peaksLR[index]["0.511 MeV"].first != -10){ + for(int i = 1; i < rangesLR[index]->size(); i++){ + TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesLR[index]->at(i)) - histo->FindBin(rangesLR[index]-> at(i-1)), rangesLR[index]-> at(i-1), rangesLR[index]->at(i)); + int j = 1; + for (int bin = histo->FindBin(rangesLR[index]->at(i-1)) ; bin < histo -> FindBin(rangesLR[index]->at(i)) +1 ; bin++){ + binHisto -> SetBinContent( j, histo->GetBinContent(bin)); + j++; + } + energyBinLR[index][i] = binHisto -> GetMean(); + binHisto -> Delete(); + } + } + } - energyBinLR[index][i] = binHisto -> GetMean(); - binHisto -> Delete(); - } - } - - } - if(!source.compare(Co60)){ - rangesLR[index] = new std::vector; - peaksLR[index] = Co60SpectrumAnalyzer(histo,rangesLR[index]); - + if(!source.compare(Co60)){ + rangesLR[index] = new std::vector; + peaksLR[index] = Co60SpectrumAnalyzer(histo,rangesLR[index]); + if (peaksLR[index]["1.173 MeV"].first > 0){ + for(auto peak : peaksLR[index] ) + { + histo -> GetXaxis() -> SetRangeUser(0.,2.*peak.second.first); + break; + } + } - if (peaksLR[index]["1.173 MeV"].first== -9999){ - peaksLR.erase(index); - rangesLR.erase(index); - peaksLR[index]["1.173 MeV"].first = -10; - peaksLR[index]["1.332 MeV"].first = -10; - - } + if (peaksLR[index]["1.173 MeV"].first== -9999){ + for(auto peak : peaksLR[index] ) + { + histo -> GetXaxis() -> SetRangeUser(0.,50); + break; + } + peaksLR.erase(index); + rangesLR.erase(index); + peaksLR[index]["1.173 MeV"].first = -10; + peaksLR[index]["1.332 MeV"].first = -10; + + } + + if( rangesLR[index] -> size() > 10){ + energy511LR = peaksLR[index]["1.173 MeV"].first; + energy1275LR = peaksLR[index]["1.332 MeV"].first; + theIndex = index; + energyBinValCo60[0] = rangesLR[index] -> at(7); + energyBinValCo60[1] = rangesLR[index] -> at(8); + energyBinValCo60[2] = rangesLR[index]-> at(9); + outTrees[index] -> Fill(); + } + + + if(peaksLR[index]["1.173 MeV"].first != -10){ + + for(int i = 1; i < rangesLR[index]->size(); i++){ + TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesLR[index]->at(i)) - histo->FindBin(rangesLR[index]-> at(i-1)), rangesLR[index]-> at(i-1), rangesLR[index]->at(i)); + int j = 1; + for (int bin = histo->FindBin(rangesLR[index]->at(i-1)) ; bin < histo -> FindBin(rangesLR[index]->at(i)) +1 ; bin++){ + binHisto -> SetBinContent( j, histo->GetBinContent(bin)); + j++; + } + energyBinLR[index][i] = binHisto -> GetMean(); + + binHisto -> Delete(); + } + } + + } - energy511LR = peaksLR[index]["1.173 MeV"].first; - energy1275LR = peaksLR[index]["1.332 MeV"].first; - theIndex = index; - outTrees[index] -> Fill(); - - - - for(auto peak : peaksLR[index] ) - { - histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first); - break; - } - - - histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum()); - if(peaksLR[index]["1.173 MeV"].first != -10){ - for(int i = 1; i < rangesLR[index]->size(); i++){ - TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesLR[index]->at(i)) - histo->FindBin(rangesLR[index]-> at(i-1)), rangesLR[index]-> at(i-1), rangesLR[index]->at(i)); - int j = 1; - for (int bin = histo->FindBin(rangesLR[index]->at(i-1)) ; bin < histo -> FindBin(rangesLR[index]->at(i)) +1 ; bin++){ - binHisto -> SetBinContent( j, histo->GetBinContent(bin)); - j++; - } - energyBinLR[index][i] = binHisto -> GetMean(); - binHisto -> Delete(); + latex -> Draw("same"); + histo -> Write(); + c -> Print(Form("%s/energy/c_energy__%s.png",plotDir.c_str(),label.c_str())); + c -> Print(Form("%s/energy/c_energy__%s.pdf",plotDir.c_str(),label.c_str())); + delete c; } - } - - } - - if( source.compare(Na22) && source.compare(Co60)){ - std::cout << " Missspelled radioactive source " << std::endl; - return(0); + } } - - - - - - - latex -> Draw("same"); - histo -> Write(); - c -> Print(Form("%s/energy/c_energy__%s.png",plotDir.c_str(),label.c_str())); - c -> Print(Form("%s/energy/c_energy__%s.pdf",plotDir.c_str(),label.c_str())); - delete c; - - } - } - } @@ -641,7 +651,7 @@ int main(int argc, char** argv) - + //------------------------ @@ -921,7 +931,7 @@ int main(int argc, char** argv) energyRatioCorr = fitFunc_energyRatioCorr[index2]->Eval(anEvent->energyR/anEvent->energyL) - fitFunc_energyRatioCorr[index2]->Eval(fitFunc_energyRatio[index2]->GetParameter(1)); - + if( h1_deltaT_energyRatioCorr[index2] == NULL ) { @@ -963,10 +973,10 @@ int main(int argc, char** argv) { double index2( 10000000*iEnergyBin + index1 ); - outTrees2[index2] = new TTree(Form("dataRes_bar%02dL-R_Vov%.1f_th%02.0f_enBin%02d",iBar,Vov,vth1,iEnergyBin),Form("dataRes_bar%02dL-R_Vov%.1f_th%02.0f_enBin%02d",iBar,Vov,vth1,iEnergyBin)); - outTrees2[index2] -> Branch("energyBin", &enBin); - outTrees2[index2] -> Branch("timeResolution",&timeRes); - outTrees2[index2] -> Branch("indexID2",&theIndex2); + outTrees2[index2] = new TTree(Form("dataRes_bar%02dL-R_Vov%.1f_th%02.0f_enBin%02d",iBar,Vov,vth1,iEnergyBin),Form("dataRes_bar%02dL-R_Vov%.1f_th%02.0f_enBin%02d",iBar,Vov,vth1,iEnergyBin)); + outTrees2[index2] -> Branch("energyBin", &enBin); + outTrees2[index2] -> Branch("timeResolution",&timeRes); + outTrees2[index2] -> Branch("indexID2",&theIndex2); diff --git a/main/moduleCharacterization_step3.cpp b/main/moduleCharacterization_step3.cpp index e4063cce..2c28f84c 100644 --- a/main/moduleCharacterization_step3.cpp +++ b/main/moduleCharacterization_step3.cpp @@ -51,8 +51,8 @@ int main(int argc, char** argv) typedef std::numeric_limits dbl; std::cout.precision(dbl::max_digits10); - -//--- get parameters + + //--- get parameters std::string plotDir = opts.GetOpt("Output.plotDirStep3"); //system(Form("rm -r %s", plotDir.c_str())); system(Form("mkdir -p %s",plotDir.c_str())); @@ -71,6 +71,7 @@ int main(int argc, char** argv) //--- open files and make the tree chain //--- Prima di eseguire cambia il numero del run nel config file + std::string inputDir1 = opts.GetOpt("Input.inputDir1"); std::string inputDir = opts.GetOpt("Input.inputDir"); std::string fileBaseName1 = opts.GetOpt("Input.fileBaseName1"); std::string fileBaseName2 = opts.GetOpt("Input.fileBaseName2"); @@ -100,7 +101,7 @@ int main(int argc, char** argv) { std::string fileNameStep1; std::string fileNameStep2; - fileNameStep1 = Form("%s%s1_%s%04d.root",inputDir.c_str(),fileBaseName1.c_str(),fileBaseName2.c_str(),run); + fileNameStep1 = Form("%s%s1_%s%04d.root",inputDir1.c_str(),fileBaseName1.c_str(),fileBaseName2.c_str(),run); fileNameStep2 = Form("%s%s2_%s%04d.root",inputDir.c_str(),fileBaseName1.c_str(),fileBaseName2.c_str(),run); std::cout << ">>> Adding file " << fileNameStep1 << std::endl; std::cout << ">>> Adding file " << fileNameStep2 << std::endl; @@ -113,6 +114,7 @@ int main(int argc, char** argv) //--- define output root file std::string outFileName = opts.GetOpt("Output.outFileNameStep3"); TFile* outFile = TFile::Open(outFileName.c_str(),"RECREATE"); + //--- get plot settings TCanvas* c; @@ -123,7 +125,9 @@ int main(int argc, char** argv) std::vector LRLabels; LRLabels.push_back("L"); LRLabels.push_back("R"); - + + float refTh = opts.GetOpt("Plots.refTh"); + std::map trees; std::map trees2; std::map VovLabels; @@ -133,11 +137,13 @@ int main(int argc, char** argv) std::map map_Vovs; std::map map_ths; std::map map_EnBin; + + std::vector vec_th1; //define TGraph std::map g_tot_vs_th; std::map g_tot_vs_Vov; - std::map g_tot_vs_iBar; + std::map g_tot_vs_bar; for(auto file: listStep1){ @@ -169,11 +175,17 @@ int main(int argc, char** argv) std::string string_th = tokens[4]; string_th.erase(0,2); map_ths[stepLabel] = atof(string_th.c_str()); + float vth1 = map_ths[stepLabel]; + vec_th1.push_back(vth1); } } std::sort(stepLabels.begin(),stepLabels.end()); stepLabels.erase(std::unique(stepLabels.begin(),stepLabels.end()),stepLabels.end()); - + std::sort(vec_th1.begin(),vec_th1.end()); + vec_th1.erase(std::unique(vec_th1.begin(),vec_th1.end()),vec_th1.end()); + + + for(auto stepLabel : stepLabels) @@ -193,6 +205,7 @@ int main(int argc, char** argv) for(auto LRLabel : LRLabels ) { + //label histo std::string label(Form("bar%02d%s_%s",iBar,LRLabel.c_str(),stepLabel.c_str())); @@ -222,11 +235,13 @@ int main(int argc, char** argv) histo -> Fit(fitFunc1,"QNS+", " ", fitFunc1->GetParameter(1)-fitFunc1->GetParameter(2), fitFunc1->GetParameter(1)+fitFunc1->GetParameter(2)); fitFunc1 -> SetLineColor(kBlack); fitFunc1 -> SetLineWidth(3); - fitFunc1 -> Draw("same"); + fitFunc1 -> Draw("same"); + //outFile -> cd(); //histo -> Write(); - c -> Print(Form("%s/summaryPlots/tot/c_tot__%s.png",plotDir.c_str(),label.c_str())); + + /*c -> Print(Form("%s/summaryPlots/tot/c_tot__%s.png",plotDir.c_str(),label.c_str())); c -> Print(Form("%s/summaryPlots/tot/c_tot__%s.pdf",plotDir.c_str(),label.c_str())); - delete c; + delete c;*/ if(LRLabel=="L"){ @@ -247,13 +262,13 @@ int main(int argc, char** argv) g_tot_vs_Vov[Form("bar%02dL_%s",iBar,thLabel.c_str())] -> SetPoint(g_tot_vs_Vov[Form("bar%02dL_%s",iBar,thLabel.c_str())]->GetN(),Vov,fitFunc1->GetMaximumX()); g_tot_vs_Vov[Form("bar%02dL_%s",iBar,thLabel.c_str())] -> SetPointError(g_tot_vs_Vov[Form("bar%02dL_%s",iBar,thLabel.c_str())]->GetN()-1,0,0); - //g_tot_vs_iBar + //g_tot_vs_bar - if( g_tot_vs_iBar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] == NULL ){ - g_tot_vs_iBar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] = new TGraphErrors(); + if( g_tot_vs_bar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] == NULL ){ + g_tot_vs_bar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] = new TGraphErrors(); } - g_tot_vs_iBar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPoint(g_tot_vs_iBar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN(),iBar,fitFunc1->GetMaximumX()); - g_tot_vs_iBar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPointError(g_tot_vs_iBar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN()-1,0,0); + g_tot_vs_bar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPoint(g_tot_vs_bar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN(),iBar,fitFunc1->GetMaximumX()); + g_tot_vs_bar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPointError(g_tot_vs_bar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN()-1,0,0); } if(LRLabel=="R"){ @@ -274,19 +289,20 @@ int main(int argc, char** argv) g_tot_vs_Vov[Form("bar%02dR_%s",iBar,thLabel.c_str())] -> SetPoint(g_tot_vs_Vov[Form("bar%02dR_%s",iBar,thLabel.c_str())]->GetN(),Vov,fitFunc1->GetMaximumX()); g_tot_vs_Vov[Form("bar%02dR_%s",iBar,thLabel.c_str())] -> SetPointError(g_tot_vs_Vov[Form("bar%02dR_%s",iBar,thLabel.c_str())]->GetN()-1,0,0); - //g_tot_vs_iBar + //g_tot_vs_bar - if( g_tot_vs_iBar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] == NULL ){ - g_tot_vs_iBar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] = new TGraphErrors(); + if( g_tot_vs_bar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] == NULL ){ + g_tot_vs_bar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] = new TGraphErrors(); } - g_tot_vs_iBar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPoint(g_tot_vs_iBar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN(),iBar,fitFunc1->GetMaximumX()); - g_tot_vs_iBar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPointError(g_tot_vs_iBar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN()-1,0,0); + g_tot_vs_bar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPoint(g_tot_vs_bar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN(),iBar,fitFunc1->GetMaximumX()); + g_tot_vs_bar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPointError(g_tot_vs_bar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN()-1,0,0); } } } } - } - + } + + //------------------ //--- draw summary plots @@ -322,11 +338,15 @@ int main(int argc, char** argv) g_totL -> SetMarkerColor(1+iter); g_totL -> SetMarkerStyle(20); g_totL -> Draw("PL,same"); + outFile -> cd(); + g_totL -> Write(Form("g_tot_vs_th_bar%02dL_%s",iBar,mapIt.first.c_str())); g_totR -> SetLineColor(1+iter); g_totR -> SetMarkerColor(1+iter); g_totR -> SetMarkerStyle(25); g_totR -> Draw("PL,same"); + outFile -> cd(); + g_totR -> Write(Form("g_tot_vs_th_bar%02dR_%s",iBar,mapIt.first.c_str())); latex = new TLatex(0.55,0.85-0.04*iter,Form("%s",mapIt.first.c_str())); latex -> SetNDC(); @@ -374,11 +394,15 @@ int main(int argc, char** argv) g_totL -> SetMarkerColor(1+iter); g_totL -> SetMarkerStyle(20); g_totL -> Draw("PL,same"); + outFile -> cd(); + g_totL -> Write(Form("g_tot_vs_Vov_bar%02dL_%s",iBar,mapIt.first.c_str())); g_totR -> SetLineColor(1+iter); g_totR -> SetMarkerColor(1+iter); g_totR -> SetMarkerStyle(25); g_totR -> Draw("PL,same"); + outFile -> cd(); + g_totR -> Write(Form("g_tot_vs_Vov_bar%02dR_%s",iBar,mapIt.first.c_str())); latex = new TLatex(0.55,0.85-0.04*iter,Form("%s",mapIt.first.c_str())); latex -> SetNDC(); @@ -398,6 +422,10 @@ int main(int argc, char** argv) for(auto mapIt1 : thLabels){ + //std::cout<<"ref"<GetN()>0){ g_totL -> SetLineColor(kBlue); g_totL -> SetMarkerColor(kBlue); g_totL -> SetMarkerStyle(20); g_totL -> Draw("P"); + outFile -> cd(); + g_totL -> Write(Form("g_tot_vs_barL_%s_%s",mapIt2.first.c_str(),mapIt1.first.c_str())); - std::string label_prova = "Blue: L Red: R"; + /*std::string label_prova = "Blue: L Red: R"; latex = new TLatex(0.55,0.85-0.04,label_prova.c_str()); latex -> SetNDC(); latex -> SetTextFont(42); latex -> SetTextSize(0.04); latex -> SetTextColor(kBlack); - latex -> Draw("same"); + latex -> Draw("same");*/ + TLegend* legenda = new TLegend(0.55,0.80,0.70,0.85); + legenda->AddEntry(g_totL, "Left", "P"); + legenda -> SetBorderSize(0); + legenda->Draw(); + c1 -> Print(Form("%s/summaryPlots/tot/c_tot_vs_bar_%s_%s.png",plotDir.c_str(),mapIt2.first.c_str(),mapIt1.first.c_str())); c1 -> Print(Form("%s/summaryPlots/tot/c_tot_vs_bar_%s_%s.pdf",plotDir.c_str(),mapIt2.first.c_str(),mapIt1.first.c_str())); } } - if(g_tot_vs_iBar[Form("R_%s_%s",mapIt2.first.c_str(),mapIt1.first.c_str())]){ - TGraphErrors* g_totR = g_tot_vs_iBar[label2]; + if(g_tot_vs_bar[Form("R_%s_%s",mapIt2.first.c_str(),mapIt1.first.c_str())]){ + TGraphErrors* g_totR = g_tot_vs_bar[label2]; if(g_totR->GetN()>0){ g_totR -> SetLineColor(kRed); g_totR -> SetMarkerColor(kRed); g_totR -> SetMarkerStyle(20); g_totR -> Draw("P,same"); + outFile -> cd(); + g_totR -> Write(Form("g_tot_vs_barR_%s_%s",mapIt2.first.c_str(),mapIt1.first.c_str())); - + TLegend* legenda = new TLegend(0.55,0.85,0.70,0.89); + legenda->AddEntry(g_totR, "Right", "P"); + legenda -> SetBorderSize(0); + legenda->Draw(); c1 -> Print(Form("%s/summaryPlots/tot/c_tot_vs_bar_%s_%s.png",plotDir.c_str(),mapIt2.first.c_str(),mapIt1.first.c_str())); @@ -454,7 +494,6 @@ int main(int argc, char** argv) - //int d=0; std::vector vec_Vov; std::vector vec_th; @@ -545,31 +584,31 @@ int main(int argc, char** argv) std::map g_en511_vs_th; std::map g_en511_vs_Vov; - std::map g_en511_vs_iBar; + std::map g_en511_vs_bar; std::map g_en1275_vs_th; std::map g_en1275_vs_Vov; - std::map g_en1275_vs_iBar; + std::map g_en1275_vs_bar; std::map g_enRatio_vs_th; std::map g_enRatio_vs_Vov; - std::map g_enRatio_vs_iBar; + std::map g_enRatio_vs_bar; std::map g_en511_vs_th_L; std::map g_en511_vs_Vov_L; - std::map g_en511_vs_iBar_L; + std::map g_en511_vs_bar_L; std::map g_en1275_vs_th_L; std::map g_en1275_vs_Vov_L; - std::map g_en1275_vs_iBar_L; + std::map g_en1275_vs_bar_L; std::map g_enRatio_vs_th_L; std::map g_enRatio_vs_Vov_L; - std::map g_enRatio_vs_iBar_L; + std::map g_enRatio_vs_bar_L; std::map g_en511_vs_th_R; std::map g_en511_vs_Vov_R; - std::map g_en511_vs_iBar_R; + std::map g_en511_vs_bar_R; std::map g_en1275_vs_th_R; std::map g_en1275_vs_Vov_R; - std::map g_en1275_vs_iBar_R; + std::map g_en1275_vs_bar_R; std::map g_enRatio_vs_th_R; std::map g_enRatio_vs_Vov_R; - std::map g_enRatio_vs_iBar_R; + std::map g_enRatio_vs_bar_R; //int index( (10000*int(Vov*100.)) + (100*vth1) + iBar ); @@ -650,28 +689,16 @@ int main(int argc, char** argv) if( g_en511_vs_Vov[th_iBar_ID] == NULL ){ g_en511_vs_Vov[th_iBar_ID] = new TGraphErrors(); } - if( g_en511_vs_iBar[Vov_th_ID] == NULL ){ - g_en511_vs_iBar[Vov_th_ID] = new TGraphErrors(); - } - - /*if( g_en511_vs_th_L[Vov_iBar_ID] == NULL ){ - g_en511_vs_th_L[Vov_iBar_ID] = new TGraphErrors(); - } - if( g_en511_vs_Vov_L[th_iBar_ID] == NULL ){ - g_en511_vs_Vov_L[th_iBar_ID] = new TGraphErrors(); - }*/ - if( g_en511_vs_iBar_L[Vov_th_ID] == NULL ){ - g_en511_vs_iBar_L[Vov_th_ID] = new TGraphErrors(); + if( g_en511_vs_bar[Vov_th_ID] == NULL ){ + g_en511_vs_bar[Vov_th_ID] = new TGraphErrors(); } - /*if( g_en511_vs_th_R[Vov_iBar_ID] == NULL ){ - g_en511_vs_th_R[Vov_iBar_ID] = new TGraphErrors(); + if( g_en511_vs_bar_L[Vov_th_ID] == NULL ){ + g_en511_vs_bar_L[Vov_th_ID] = new TGraphErrors(); } - if( g_en511_vs_Vov_R[th_iBar_ID] == NULL ){ - g_en511_vs_Vov_R[th_iBar_ID] = new TGraphErrors(); - }*/ - if( g_en511_vs_iBar_R[Vov_th_ID] == NULL ){ - g_en511_vs_iBar_R[Vov_th_ID] = new TGraphErrors(); + + if( g_en511_vs_bar_R[Vov_th_ID] == NULL ){ + g_en511_vs_bar_R[Vov_th_ID] = new TGraphErrors(); } if (index->second>0.1){ @@ -684,9 +711,9 @@ int main(int argc, char** argv) //std::cout<<"N"<GetN()<<"Vov"<second<<"thIbarID"<SetPointError(g_en511_vs_Vov[th_iBar_ID]->GetN()-1, 0., 0.); - g_en511_vs_iBar[Vov_th_ID]->SetPoint(g_en511_vs_iBar[Vov_th_ID]->GetN(), iBar, index->second); - //std::cout<<"N"<GetN()<<"Vov"<second<<"thVovID"<SetPointError(g_en511_vs_iBar[Vov_th_ID]->GetN()-1, 0., 0.); + g_en511_vs_bar[Vov_th_ID]->SetPoint(g_en511_vs_bar[Vov_th_ID]->GetN(), iBar, index->second); + //std::cout<<"N"<GetN()<<"Vov"<second<<"thVovID"<SetPointError(g_en511_vs_bar[Vov_th_ID]->GetN()-1, 0., 0.); } @@ -698,8 +725,8 @@ int main(int argc, char** argv) g_en511_vs_Vov_L[th_iBar_ID]->SetPoint(g_en511_vs_Vov_L[th_iBar_ID]->GetN(), Vov, map_en511L[index->first]); g_en511_vs_Vov_L[th_iBar_ID]->SetPointError(g_en511_vs_Vov_L[th_iBar_ID]->GetN()-1, 0., 0.);*/ - g_en511_vs_iBar_L[Vov_th_ID]->SetPoint(g_en511_vs_iBar_L[Vov_th_ID]->GetN(), iBar, map_en511L[index->first]); - g_en511_vs_iBar_L[Vov_th_ID]->SetPointError(g_en511_vs_iBar_L[Vov_th_ID]->GetN()-1, 0., 0.); + g_en511_vs_bar_L[Vov_th_ID]->SetPoint(g_en511_vs_bar_L[Vov_th_ID]->GetN(), iBar, map_en511L[index->first]); + g_en511_vs_bar_L[Vov_th_ID]->SetPointError(g_en511_vs_bar_L[Vov_th_ID]->GetN()-1, 0., 0.); } @@ -711,8 +738,8 @@ int main(int argc, char** argv) g_en511_vs_Vov_R[th_iBar_ID]->SetPoint(g_en511_vs_Vov_R[th_iBar_ID]->GetN(), Vov, map_en511R[index->first]); g_en511_vs_Vov_R[th_iBar_ID]->SetPointError(g_en511_vs_Vov_R[th_iBar_ID]->GetN()-1, 0., 0.);*/ - g_en511_vs_iBar_R[Vov_th_ID]->SetPoint(g_en511_vs_iBar_R[Vov_th_ID]->GetN(), iBar, map_en511R[index->first]); - g_en511_vs_iBar_R[Vov_th_ID]->SetPointError(g_en511_vs_iBar_R[Vov_th_ID]->GetN()-1, 0., 0.); + g_en511_vs_bar_R[Vov_th_ID]->SetPoint(g_en511_vs_bar_R[Vov_th_ID]->GetN(), iBar, map_en511R[index->first]); + g_en511_vs_bar_R[Vov_th_ID]->SetPointError(g_en511_vs_bar_R[Vov_th_ID]->GetN()-1, 0., 0.); } @@ -741,26 +768,16 @@ int main(int argc, char** argv) if( g_en1275_vs_Vov[th_iBar_ID] == NULL ){ g_en1275_vs_Vov[th_iBar_ID] = new TGraphErrors(); } - if( g_en1275_vs_iBar[Vov_th_ID] == NULL ){ - g_en1275_vs_iBar[Vov_th_ID] = new TGraphErrors(); - } - /*if( g_en1275_vs_th_L[Vov_iBar_ID] == NULL ){ - g_en1275_vs_th_L[Vov_iBar_ID] = new TGraphErrors(); + if( g_en1275_vs_bar[Vov_th_ID] == NULL ){ + g_en1275_vs_bar[Vov_th_ID] = new TGraphErrors(); } - if( g_en1275_vs_Vov_L[th_iBar_ID] == NULL ){ - g_en1275_vs_Vov_L[th_iBar_ID] = new TGraphErrors(); - }*/ - if( g_en1275_vs_iBar_L[Vov_th_ID] == NULL ){ - g_en1275_vs_iBar_L[Vov_th_ID] = new TGraphErrors(); - } - /*if( g_en1275_vs_th_R[Vov_iBar_ID] == NULL ){ - g_en1275_vs_th_R[Vov_iBar_ID] = new TGraphErrors(); + + if( g_en1275_vs_bar_L[Vov_th_ID] == NULL ){ + g_en1275_vs_bar_L[Vov_th_ID] = new TGraphErrors(); } - if( g_en1275_vs_Vov_R[th_iBar_ID] == NULL ){ - g_en1275_vs_Vov_R[th_iBar_ID] = new TGraphErrors(); - }*/ - if( g_en1275_vs_iBar_R[Vov_th_ID] == NULL ){ - g_en1275_vs_iBar_R[Vov_th_ID] = new TGraphErrors(); + + if( g_en1275_vs_bar_R[Vov_th_ID] == NULL ){ + g_en1275_vs_bar_R[Vov_th_ID] = new TGraphErrors(); } if (index->second>0.1){ @@ -771,8 +788,8 @@ int main(int argc, char** argv) g_en1275_vs_Vov[th_iBar_ID]->SetPoint(g_en1275_vs_Vov[th_iBar_ID]->GetN(), Vov, index->second); g_en1275_vs_Vov[th_iBar_ID]->SetPointError(g_en1275_vs_Vov[th_iBar_ID]->GetN()-1, 0., 0.); - g_en1275_vs_iBar[Vov_th_ID]->SetPoint(g_en1275_vs_iBar[Vov_th_ID]->GetN(), iBar, index->second); - g_en1275_vs_iBar[Vov_th_ID]->SetPointError(g_en1275_vs_iBar[Vov_th_ID]->GetN()-1, 0., 0.); + g_en1275_vs_bar[Vov_th_ID]->SetPoint(g_en1275_vs_bar[Vov_th_ID]->GetN(), iBar, index->second); + g_en1275_vs_bar[Vov_th_ID]->SetPointError(g_en1275_vs_bar[Vov_th_ID]->GetN()-1, 0., 0.); } if (map_en1275R[index->first]>0.1){ @@ -783,8 +800,8 @@ int main(int argc, char** argv) g_en1275_vs_Vov_R[th_iBar_ID]->SetPoint(g_en1275_vs_Vov_R[th_iBar_ID]->GetN(), Vov, map_en1275R[index->first]); g_en1275_vs_Vov_R[th_iBar_ID]->SetPointError(g_en1275_vs_Vov_R[th_iBar_ID]->GetN()-1, 0., 0.);*/ - g_en1275_vs_iBar_R[Vov_th_ID]->SetPoint(g_en1275_vs_iBar_R[Vov_th_ID]->GetN(), iBar, map_en1275R[index->first]); - g_en1275_vs_iBar_R[Vov_th_ID]->SetPointError(g_en1275_vs_iBar_R[Vov_th_ID]->GetN()-1, 0., 0.); + g_en1275_vs_bar_R[Vov_th_ID]->SetPoint(g_en1275_vs_bar_R[Vov_th_ID]->GetN(), iBar, map_en1275R[index->first]); + g_en1275_vs_bar_R[Vov_th_ID]->SetPointError(g_en1275_vs_bar_R[Vov_th_ID]->GetN()-1, 0., 0.); } if (map_en1275L[index->first]>0.1){ @@ -795,8 +812,8 @@ int main(int argc, char** argv) g_en1275_vs_Vov_L[th_iBar_ID]->SetPoint(g_en1275_vs_Vov_L[th_iBar_ID]->GetN(), Vov, map_en1275L[index->first]); g_en1275_vs_Vov_L[th_iBar_ID]->SetPointError(g_en1275_vs_Vov_L[th_iBar_ID]->GetN()-1, 0., 0.);*/ - g_en1275_vs_iBar_L[Vov_th_ID]->SetPoint(g_en1275_vs_iBar_L[Vov_th_ID]->GetN(), iBar, map_en1275L[index->first]); - g_en1275_vs_iBar_L[Vov_th_ID]->SetPointError(g_en1275_vs_iBar_L[Vov_th_ID]->GetN()-1, 0., 0.); + g_en1275_vs_bar_L[Vov_th_ID]->SetPoint(g_en1275_vs_bar_L[Vov_th_ID]->GetN(), iBar, map_en1275L[index->first]); + g_en1275_vs_bar_L[Vov_th_ID]->SetPointError(g_en1275_vs_bar_L[Vov_th_ID]->GetN()-1, 0., 0.); } } @@ -824,30 +841,16 @@ int main(int argc, char** argv) g_enRatio_vs_Vov[th_iBar_ID] = new TGraphErrors(); } - if( g_enRatio_vs_iBar[Vov_th_ID] == NULL ){ - g_enRatio_vs_iBar[Vov_th_ID] = new TGraphErrors(); + if( g_enRatio_vs_bar[Vov_th_ID] == NULL ){ + g_enRatio_vs_bar[Vov_th_ID] = new TGraphErrors(); } - /*if( g_enRatio_vs_th_L[Vov_iBar_ID] == NULL ){ - g_enRatio_vs_th_L[Vov_iBar_ID] = new TGraphErrors(); + if( g_enRatio_vs_bar_L[Vov_th_ID] == NULL ){ + g_enRatio_vs_bar_L[Vov_th_ID] = new TGraphErrors(); } - if( g_enRatio_vs_Vov_L[th_iBar_ID] == NULL ){ - g_enRatio_vs_Vov_L[th_iBar_ID] = new TGraphErrors(); - }*/ - - if( g_enRatio_vs_iBar_L[Vov_th_ID] == NULL ){ - g_enRatio_vs_iBar_L[Vov_th_ID] = new TGraphErrors(); - } - - /* if( g_enRatio_vs_th_R[Vov_iBar_ID] == NULL ){ - g_enRatio_vs_th_R[Vov_iBar_ID] = new TGraphErrors(); - } - if( g_enRatio_vs_Vov_R[th_iBar_ID] == NULL ){ - g_enRatio_vs_Vov_R[th_iBar_ID] = new TGraphErrors(); - }*/ - - if( g_enRatio_vs_iBar_R[Vov_th_ID] == NULL ){ - g_enRatio_vs_iBar_R[Vov_th_ID] = new TGraphErrors(); + + if( g_enRatio_vs_bar_R[Vov_th_ID] == NULL ){ + g_enRatio_vs_bar_R[Vov_th_ID] = new TGraphErrors(); } if (index->second>0.1){ @@ -858,8 +861,8 @@ int main(int argc, char** argv) g_enRatio_vs_Vov[th_iBar_ID]->SetPoint(g_enRatio_vs_Vov[th_iBar_ID]->GetN(), Vov, index->second); g_enRatio_vs_Vov[th_iBar_ID]->SetPointError(g_enRatio_vs_Vov[th_iBar_ID]->GetN()-1, 0., 0.); - g_enRatio_vs_iBar[Vov_th_ID]->SetPoint(g_enRatio_vs_iBar[Vov_th_ID]->GetN(), iBar, index->second); - g_enRatio_vs_iBar[Vov_th_ID]->SetPointError(g_enRatio_vs_iBar[Vov_th_ID]->GetN()-1, 0., 0.); + g_enRatio_vs_bar[Vov_th_ID]->SetPoint(g_enRatio_vs_bar[Vov_th_ID]->GetN(), iBar, index->second); + g_enRatio_vs_bar[Vov_th_ID]->SetPointError(g_enRatio_vs_bar[Vov_th_ID]->GetN()-1, 0., 0.); } if (map_enRatioR[index->first]>0.1){ @@ -870,8 +873,8 @@ int main(int argc, char** argv) g_enRatio_vs_Vov_R[th_iBar_ID]->SetPoint(g_enRatio_vs_Vov_R[th_iBar_ID]->GetN(), Vov, map_enRatioR[index->first]); g_enRatio_vs_Vov_R[th_iBar_ID]->SetPointError(g_enRatio_vs_Vov_R[th_iBar_ID]->GetN()-1, 0., 0.);*/ - g_enRatio_vs_iBar_R[Vov_th_ID]->SetPoint(g_enRatio_vs_iBar_R[Vov_th_ID]->GetN(), iBar, map_enRatioR[index->first]); - g_enRatio_vs_iBar_R[Vov_th_ID]->SetPointError(g_enRatio_vs_iBar_R[Vov_th_ID]->GetN()-1, 0., 0.); + g_enRatio_vs_bar_R[Vov_th_ID]->SetPoint(g_enRatio_vs_bar_R[Vov_th_ID]->GetN(), iBar, map_enRatioR[index->first]); + g_enRatio_vs_bar_R[Vov_th_ID]->SetPointError(g_enRatio_vs_bar_R[Vov_th_ID]->GetN()-1, 0., 0.); } if (map_enRatioL[index->first]>0.1){ @@ -882,44 +885,43 @@ int main(int argc, char** argv) g_enRatio_vs_Vov_L[th_iBar_ID]->SetPoint(g_enRatio_vs_Vov_L[th_iBar_ID]->GetN(), Vov, map_enRatioL[index->first]); g_enRatio_vs_Vov_L[th_iBar_ID]->SetPointError(g_enRatio_vs_Vov_L[th_iBar_ID]->GetN()-1, 0., 0.);*/ - g_enRatio_vs_iBar_L[Vov_th_ID]->SetPoint(g_enRatio_vs_iBar_L[Vov_th_ID]->GetN(), iBar, map_enRatioL[index->first]); - g_enRatio_vs_iBar_L[Vov_th_ID]->SetPointError(g_enRatio_vs_iBar_L[Vov_th_ID]->GetN()-1, 0., 0.); + g_enRatio_vs_bar_L[Vov_th_ID]->SetPoint(g_enRatio_vs_bar_L[Vov_th_ID]->GetN(), iBar, map_enRatioL[index->first]); + g_enRatio_vs_bar_L[Vov_th_ID]->SetPointError(g_enRatio_vs_bar_L[Vov_th_ID]->GetN()-1, 0., 0.); } } - - - - - + + + + std::map c_en511_vs_th; std::map c_en511_vs_Vov; - std::map c_en511_vs_iBar; + std::map c_en511_vs_bar; std::map c_en1275_vs_th; std::map c_en1275_vs_Vov; - std::map c_en1275_vs_iBar; + std::map c_en1275_vs_bar; std::map c_enRatio_vs_th; std::map c_enRatio_vs_Vov; - std::map c_enRatio_vs_iBar; + std::map c_enRatio_vs_bar; std::map c_en511_vs_th_L; std::map c_en511_vs_Vov_L; - std::map c_en511_vs_iBar_L; + std::map c_en511_vs_bar_L; std::map c_en1275_vs_th_L; std::map c_en1275_vs_Vov_L; - std::map c_en1275_vs_iBar_L; + std::map c_en1275_vs_bar_L; std::map c_enRatio_vs_th_L; std::map c_enRatio_vs_Vov_L; - std::map c_enRatio_vs_iBar_L; + std::map c_enRatio_vs_bar_L; std::map c_en511_vs_th_R; std::map c_en511_vs_Vov_R; - std::map c_en511_vs_iBar_R; + std::map c_en511_vs_bar_R; std::map c_en1275_vs_th_R; std::map c_en1275_vs_Vov_R; - std::map c_en1275_vs_iBar_R; + std::map c_en1275_vs_bar_R; std::map c_enRatio_vs_th_R; std::map c_enRatio_vs_Vov_R; - std::map c_enRatio_vs_iBar_R; + std::map c_enRatio_vs_bar_R; std::map iter; @@ -1018,6 +1020,8 @@ int main(int argc, char** argv) g_energy -> SetMarkerColor(1+iter[bar]); g_energy -> SetMarkerStyle(20); g_energy -> Draw("PL,same"); + outFile -> cd(); + g_energy -> Write(Form("g_en511_vs_th_Vov%.01f_bar%02d",Vov,bar)); std::string VovLabel = Form("Vov%.01f", Vov); latex = new TLatex(0.55,0.85-0.04*iter[bar],VovLabel.c_str()); @@ -1055,6 +1059,8 @@ int main(int argc, char** argv) g_energy1275 -> SetMarkerColor(1+iter[bar]); g_energy1275 -> SetMarkerStyle(20); g_energy1275 -> Draw("PL,same"); + outFile -> cd(); + g_energy1275 -> Write(Form("g_en1275_vs_th_Vov%.01f_bar%02d",Vov,bar)); latex -> Draw("same"); @@ -1085,6 +1091,9 @@ int main(int argc, char** argv) g_energyRatio -> SetMarkerColor(1+iter[bar]); g_energyRatio -> SetMarkerStyle(20); g_energyRatio -> Draw("PL,same"); + outFile -> cd(); + g_energyRatio -> Write(Form("g_enRatio_vs_th_Vov%.01f_bar%02d",Vov,bar)); + latex -> Draw("same"); @@ -1255,6 +1264,8 @@ int main(int argc, char** argv) g_energy -> SetMarkerColor(1+iter[bar]); g_energy -> SetMarkerStyle(20); g_energy -> Draw("PL,same"); + outFile -> cd(); + g_energy -> Write(Form("g_en511_vs_Vov_th%d_bar%02d",th,bar)); std::string thLabel = Form("th%d", th); latex = new TLatex(0.55,0.85-0.04*iter[bar],thLabel.c_str()); @@ -1292,6 +1303,8 @@ int main(int argc, char** argv) g_energy1275 -> SetMarkerColor(1+iter[bar]); g_energy1275 -> SetMarkerStyle(20); g_energy1275 -> Draw("PL,same"); + outFile -> cd(); + g_energy1275 -> Write(Form("g_en1275_vs_Vov_th%d_bar%02d",th,bar)); latex -> Draw("same"); @@ -1322,6 +1335,8 @@ int main(int argc, char** argv) g_energyRatio -> SetMarkerColor(1+iter[bar]); g_energyRatio -> SetMarkerStyle(20); g_energyRatio -> Draw("PL,same"); + outFile -> cd(); + g_energyRatio -> Write(Form("g_enRatio_vs_Vov_th%d_bar%02d",th,bar)); latex -> Draw("same"); @@ -1399,94 +1414,64 @@ int main(int argc, char** argv) //summary plots Energy Peak 511, 1275, Ratio vs iBar (LR, L, R) - for(std::map::iterator index = g_en511_vs_iBar.begin(); index != g_en511_vs_iBar.end(); index++){ - + for(std::map::iterator index = g_en511_vs_bar.begin(); index != g_en511_vs_bar.end(); index++){ + int th; float Vov; int Vov_th_ID; Vov = float((int(index->first/10000))/100); th = int((index->first - Vov*1000000)/100); Vov_th_ID = 10000*int(Vov*100.) + 100*th; - - if(c_en511_vs_iBar[Vov_th_ID]==NULL){ - c_en511_vs_iBar[Vov_th_ID] = new TCanvas(Form("c_en511_vs_bar_Vov%.01f_th%d",Vov,th),Form("c_en511_vs_bar_Vov%.01f_th%d",Vov,th)); + + if( th != refTh ) continue; + + if(c_en511_vs_bar[Vov_th_ID]==NULL){ + c_en511_vs_bar[Vov_th_ID] = new TCanvas(Form("c_en511_vs_bar_Vov%.01f_th%d",Vov,th),Form("c_en511_vs_bar_Vov%.01f_th%d",Vov,th)); iter[Vov_th_ID] = 0; TH1F* hPad = (TH1F*)( gPad->DrawFrame(0.,0.,17.,25.) ); hPad -> SetTitle(";ID bar;energy [a.u.]"); hPad -> Draw(); gPad -> SetGridy(); } - if(c_en511_vs_iBar_L[Vov_th_ID]==NULL){ - c_en511_vs_iBar_L[Vov_th_ID] = new TCanvas(Form("c_en511_vs_barL_Vov%.01f_th%d",Vov,th),Form("c_en511_vs_barL_Vov%.01f_th%d",Vov,th)); - TH1F* hPad = (TH1F*)( gPad->DrawFrame(0.,0.,17.,25.) ); - hPad -> SetTitle(";ID bar;energy [a.u.]"); - hPad -> Draw(); - gPad -> SetGridy(); - } - if(c_en511_vs_iBar_R[Vov_th_ID]==NULL){ - c_en511_vs_iBar_R[Vov_th_ID] = new TCanvas(Form("c_en511_vs_barR_Vov%.01f_th%d",Vov,th),Form("c_en511_vs_barR_Vov%.01f_th%d",Vov,th)); - TH1F* hPad = (TH1F*)( gPad->DrawFrame(0.,0.,17.,25.) ); - hPad -> SetTitle(";ID bar;energy [a.u.]"); - hPad -> Draw(); - gPad -> SetGridy(); - } - - if(c_en1275_vs_iBar[Vov_th_ID]==NULL){ - c_en1275_vs_iBar[Vov_th_ID] = new TCanvas(Form("c_en1275_vs_bar_Vov%.01f_th%d",Vov,th),Form("c_en1275_vs_bar_Vov%.01f_th%d",Vov,th)); - TH1F* hPad = (TH1F*)( gPad->DrawFrame(0.,0.,17.,35.) ); - hPad -> SetTitle(";ID bar;energy [a.u.]"); - hPad -> Draw(); - gPad -> SetGridy(); - } - - if(c_en1275_vs_iBar_R[Vov_th_ID]==NULL){ - c_en1275_vs_iBar_R[Vov_th_ID] = new TCanvas(Form("c_en1275_vs_barR_Vov%.01f_th%d",Vov,th),Form("c_en1275_vs_barR_Vov%.01f_th%d",Vov,th)); - TH1F* hPad = (TH1F*)( gPad->DrawFrame(0.,0.,17.,35.) ); - hPad -> SetTitle(";ID bar;energy [a.u.]"); - hPad -> Draw(); - gPad -> SetGridy(); - } - - if(c_en1275_vs_iBar_L[Vov_th_ID]==NULL){ - c_en1275_vs_iBar_L[Vov_th_ID] = new TCanvas(Form("c_en1275_vs_barL_Vov%.01f_th%d",Vov,th),Form("c_en1275_vs_barL_Vov%.01f_th%d",Vov,th)); + + if(c_en1275_vs_bar[Vov_th_ID]==NULL){ + c_en1275_vs_bar[Vov_th_ID] = new TCanvas(Form("c_en1275_vs_bar_Vov%.01f_th%d",Vov,th),Form("c_en1275_vs_bar_Vov%.01f_th%d",Vov,th)); TH1F* hPad = (TH1F*)( gPad->DrawFrame(0.,0.,17.,35.) ); hPad -> SetTitle(";ID bar;energy [a.u.]"); hPad -> Draw(); gPad -> SetGridy(); } - - if(c_enRatio_vs_iBar[Vov_th_ID]==NULL){ - c_enRatio_vs_iBar[Vov_th_ID] = new TCanvas(Form("c_enRatio_vs_bar_Vov%.01f_th%d",Vov,th),Form("c_enRatio_vs_bar_Vov%.01f_th%d",Vov,th)); - TH1F* hPad = (TH1F*)( gPad->DrawFrame(0.,1.,17.,3.) ); - hPad -> SetTitle(";ID bar;energy [a.u.]"); - hPad -> Draw(); - gPad -> SetGridy(); - } - - if(c_enRatio_vs_iBar_R[Vov_th_ID]==NULL){ - c_enRatio_vs_iBar_R[Vov_th_ID] = new TCanvas(Form("c_enRatio_vs_barR_Vov%.01f_th%d",Vov,th),Form("c_enRatio_vs_barR_Vov%.01f_th%d",Vov,th)); - TH1F* hPad = (TH1F*)( gPad->DrawFrame(0.,1.,17.,3.) ); - hPad -> SetTitle(";ID bar;energy [a.u.]"); - hPad -> Draw(); - gPad -> SetGridy(); - } - - if(c_enRatio_vs_iBar_L[Vov_th_ID]==NULL){ - c_enRatio_vs_iBar_L[Vov_th_ID] = new TCanvas(Form("c_enRatio_vs_barL_Vov%.01f_th%d",Vov,th),Form("c_enRatio_vs_barL_Vov%.01f_th%d",Vov,th)); + + if(c_enRatio_vs_bar[Vov_th_ID]==NULL){ + c_enRatio_vs_bar[Vov_th_ID] = new TCanvas(Form("c_enRatio_vs_bar_Vov%.01f_th%d",Vov,th),Form("c_enRatio_vs_bar_Vov%.01f_th%d",Vov,th)); TH1F* hPad = (TH1F*)( gPad->DrawFrame(0.,1.,17.,3.) ); - hPad -> SetTitle(";ID bar;energy [a.u.]"); + hPad -> SetTitle(";ID bar;energy_{1275 keV} / energy_{511 keV}"); hPad -> Draw(); gPad -> SetGridy(); } + + c_en511_vs_bar[Vov_th_ID]->cd(); + TGraph* g_energy = g_en511_vs_bar[Vov_th_ID]; + TGraph* g_energy_L = g_en511_vs_bar_L[Vov_th_ID]; + TGraph* g_energy_R = g_en511_vs_bar_R[Vov_th_ID]; - c_en511_vs_iBar[Vov_th_ID]->cd(); - TGraph* g_energy = g_en511_vs_iBar[Vov_th_ID]; - g_energy -> SetLineColor(1+iter[Vov_th_ID]); g_energy -> SetMarkerColor(1+iter[Vov_th_ID]); g_energy -> SetMarkerStyle(20); g_energy -> Draw("PL,same"); - + outFile -> cd(); + g_energy -> Write(Form("g_en511_vs_bar_Vov%.01f_th%d",Vov,th)); + g_energy_L -> SetLineStyle(9); + g_energy_L -> SetMarkerStyle(24); + g_energy_L -> Draw("PL,same"); + outFile -> cd(); + g_energy_L -> Write(Form("g_en511_vs_barL_Vov%.01f_th%d",Vov,th)); + g_energy_R -> SetLineStyle(9); + g_energy_R -> SetMarkerStyle(25); + g_energy_R -> Draw("PL,same"); + outFile -> cd(); + g_energy_R -> Write(Form("g_en511_vs_barR_Vov%.01f_th%d",Vov,th)); + /*std::string VovthLabel = Form("Vov%.01f th%d", Vov,th); latex = new TLatex(0.55,0.85-0.04*iter[Vov_th_ID],VovthLabel.c_str()); latex -> SetNDC(); @@ -1494,197 +1479,112 @@ int main(int argc, char** argv) latex -> SetTextSize(0.04); latex -> SetTextColor(kBlack+iter[Vov_th_ID]); latex -> Draw("same");*/ - - c_en511_vs_iBar_L[Vov_th_ID]->cd(); - TGraph* g_energy_L = g_en511_vs_iBar_L[Vov_th_ID]; - - g_energy_L -> SetLineStyle(2); - g_energy_L -> SetLineColor(1+iter[Vov_th_ID]); - g_energy_L -> SetMarkerColor(1+iter[Vov_th_ID]); - g_energy_L -> SetMarkerStyle(24); - g_energy_L -> SetMarkerSize(0.8); - g_energy_L -> Draw("PL,same"); - - //latex -> Draw("same"); - - c_en511_vs_iBar_R[Vov_th_ID]->cd(); - TGraph* g_energy_R = g_en511_vs_iBar_R[Vov_th_ID]; - - g_energy_R -> SetLineStyle(2); - g_energy_R -> SetLineColor(1+iter[Vov_th_ID]); - g_energy_R -> SetMarkerColor(1+iter[Vov_th_ID]); - g_energy_R -> SetMarkerStyle(25); - g_energy_R -> SetMarkerSize(0.8); - g_energy_R -> Draw("PL,same"); - - c_en511_vs_iBar[Vov_th_ID]->cd(); - g_energy_L -> Draw("PL,same"); - g_energy_R -> Draw("PL,same"); - - latex = new TLatex(0.20,0.87, Form("Average = %.1f, RMS = %.1f %%", g_energy ->GetMean(2), g_energy->GetRMS(2)/g_energy ->GetMean(2)*100)); + latex = new TLatex(0.25,0.87, Form("Average = %.1f, RMS = %.1f %%", g_energy ->GetMean(2), g_energy->GetRMS(2)/g_energy ->GetMean(2)*100)); latex -> SetNDC(); latex -> SetTextFont(42); latex -> SetTextSize(0.03); latex ->Draw("same"); - latex = new TLatex(0.20,0.83, Form("Left = %.1f, RMS = %.1f %%", g_energy_L ->GetMean(2), g_energy_L->GetRMS(2)/g_energy_L ->GetMean(2)*100 )); + latex = new TLatex(0.25,0.83, Form("Left = %.1f, RMS = %.1f %%", g_energy_L ->GetMean(2), g_energy_L->GetRMS(2)/g_energy_L ->GetMean(2)*100 )); latex -> SetNDC(); latex -> SetTextFont(42); latex -> SetTextSize(0.03); latex ->Draw("same"); - latex = new TLatex(0.20,0.79, Form("Right = %.1f, RMS = %.1f %%", g_energy_R ->GetMean(2), g_energy_R->GetRMS(2)/g_energy_R ->GetMean(2)*100)); + latex = new TLatex(0.25,0.79, Form("Right = %.1f, RMS = %.1f %%", g_energy_R ->GetMean(2), g_energy_R->GetRMS(2)/g_energy_R ->GetMean(2)*100)); latex -> SetNDC(); latex -> SetTextFont(42); latex -> SetTextSize(0.03); latex ->Draw("same"); - - //latex -> Draw("same"); + + TLegend* legenda = new TLegend(0.20,0.78,0.25,0.90); + legenda->AddEntry(g_energy, "", "P"); + legenda->AddEntry(g_energy_L, "", "P"); + legenda->AddEntry(g_energy_R, "", "P"); + legenda -> SetBorderSize(0); + legenda->Draw(); - c_en1275_vs_iBar[Vov_th_ID]->cd(); - TGraph* g_energy1275 = g_en1275_vs_iBar[Vov_th_ID]; + c_en1275_vs_bar[Vov_th_ID]->cd(); + TGraph* g_energy1275 = g_en1275_vs_bar[Vov_th_ID]; + TGraph* g_energy1275_L = g_en1275_vs_bar_L[Vov_th_ID]; + TGraph* g_energy1275_R = g_en1275_vs_bar_R[Vov_th_ID]; g_energy1275 -> SetLineColor(1+iter[Vov_th_ID]); g_energy1275 -> SetMarkerColor(1+iter[Vov_th_ID]); g_energy1275 -> SetMarkerStyle(20); g_energy1275 -> Draw("PL,same"); - - //latex -> Draw("same"); - - c_en1275_vs_iBar_L[Vov_th_ID]->cd(); - TGraph* g_energy1275_L = g_en1275_vs_iBar_L[Vov_th_ID]; - - g_energy1275_L -> SetLineStyle(2); - g_energy1275_L -> SetLineColor(1+iter[Vov_th_ID]); - g_energy1275_L -> SetMarkerColor(1+iter[Vov_th_ID]); - g_energy1275_L -> SetMarkerStyle(24); - g_energy1275_L -> SetMarkerSize(0.8); + outFile -> cd(); + g_energy1275 -> Write(Form("g_en1275_vs_bar_Vov%.01f_th%d",Vov,th)); + g_energy1275_L -> SetLineStyle(9); + g_energy1275_L -> SetMarkerStyle(24); g_energy1275_L -> Draw("PL,same"); - - //latex -> Draw("same"); - - c_en1275_vs_iBar_R[Vov_th_ID]->cd(); - TGraph* g_energy1275_R = g_en1275_vs_iBar_R[Vov_th_ID]; - - g_energy1275_R -> SetLineStyle(2); - g_energy1275_R -> SetLineColor(1+iter[Vov_th_ID]); - g_energy1275_R -> SetMarkerColor(1+iter[Vov_th_ID]); - g_energy1275_R -> SetMarkerStyle(25); - g_energy1275_R -> SetMarkerSize(0.8); + outFile -> cd(); + g_energy1275_L-> Write(Form("g_en1275_vs_barL_Vov%.01f_th%d",Vov,th)); + g_energy1275_R -> SetLineStyle(9); + g_energy1275_R -> SetMarkerStyle(25); g_energy1275_R -> Draw("PL,same"); - - //latex -> Draw("same"); - - c_en1275_vs_iBar[Vov_th_ID]->cd(); - g_energy1275_L -> Draw("PL,same"); - g_energy1275_R -> Draw("PL,same"); - - latex = new TLatex(0.20,0.87, Form("Average = %.1f, RMS = %.1f %%", g_energy1275 ->GetMean(2), g_energy1275->GetRMS(2)/g_energy1275 ->GetMean(2)*100)); + outFile -> cd(); + g_energy1275_R -> Write(Form("g_en1275_vs_barR_Vov%.01f_th%d",Vov,th)); + + latex = new TLatex(0.25,0.87, Form("Average = %.1f, RMS = %.1f %%", g_energy1275 ->GetMean(2), g_energy1275->GetRMS(2)/g_energy1275 ->GetMean(2)*100)); latex -> SetNDC(); latex -> SetTextFont(42); latex -> SetTextSize(0.03); latex ->Draw("same"); - latex = new TLatex(0.20,0.83, Form("Left = %.1f, RMS = %.1f %%", g_energy1275_L ->GetMean(2), g_energy1275_L->GetRMS(2)/g_energy1275_L ->GetMean(2)*100 )); + latex = new TLatex(0.25,0.83, Form("Left = %.1f, RMS = %.1f %%", g_energy1275_L ->GetMean(2), g_energy1275_L->GetRMS(2)/g_energy1275_L ->GetMean(2)*100 )); latex -> SetNDC(); latex -> SetTextFont(42); latex -> SetTextSize(0.03); latex ->Draw("same"); - latex = new TLatex(0.20,0.79, Form("Right = %.1f, RMS = %.1f %%", g_energy1275_R ->GetMean(2), g_energy1275_R->GetRMS(2)/g_energy1275_R ->GetMean(2)*100)); + latex = new TLatex(0.25,0.79, Form("Right = %.1f, RMS = %.1f %%", g_energy1275_R ->GetMean(2), g_energy1275_R->GetRMS(2)/g_energy1275_R ->GetMean(2)*100)); latex -> SetNDC(); latex -> SetTextFont(42); latex -> SetTextSize(0.03); latex ->Draw("same"); - - c_enRatio_vs_iBar[Vov_th_ID]->cd(); - TGraph* g_energyRatio = g_enRatio_vs_iBar[Vov_th_ID]; + + TLegend* legenda2 = new TLegend(0.20,0.78,0.25,0.90); + legenda2->AddEntry(g_energy1275, "", "P"); + legenda2->AddEntry(g_energy1275_L, "", "P"); + legenda2->AddEntry(g_energy1275_R, "", "P"); + legenda2 -> SetBorderSize(0); + legenda2->Draw(); + + c_enRatio_vs_bar[Vov_th_ID]->cd(); + TGraph* g_energyRatio = g_enRatio_vs_bar[Vov_th_ID]; g_energyRatio -> SetLineColor(1+iter[Vov_th_ID]); g_energyRatio -> SetMarkerColor(1+iter[Vov_th_ID]); g_energyRatio -> SetMarkerStyle(20); g_energyRatio -> Draw("PL,same"); - - //latex -> Draw("same"); - - c_enRatio_vs_iBar_L[Vov_th_ID]->cd(); - TGraph* g_energyRatio_L = g_enRatio_vs_iBar_L[Vov_th_ID]; - - g_energyRatio_L -> SetLineColor(1+iter[Vov_th_ID]); - g_energyRatio_L -> SetMarkerColor(1+iter[Vov_th_ID]); - g_energyRatio_L -> SetMarkerStyle(20); - g_energyRatio_L -> Draw("PL,same"); - - //latex -> Draw("same"); - - c_enRatio_vs_iBar_R[Vov_th_ID]->cd(); - TGraph* g_energyRatio_R = g_enRatio_vs_iBar_R[Vov_th_ID]; - - g_energyRatio_R -> SetLineColor(1+iter[Vov_th_ID]); - g_energyRatio_R -> SetMarkerColor(1+iter[Vov_th_ID]); - g_energyRatio_R -> SetMarkerStyle(20); - g_energyRatio_R -> Draw("PL,same"); - - //latex -> Draw("same"); - - + outFile -> cd(); + g_energyRatio -> Write(Form("g_enRatio_vs_bar_Vov%.01f_th%d",Vov,th)); + + ++iter[Vov_th_ID]; - } - - for(std::map::iterator index = c_en511_vs_iBar.begin(); index != c_en511_vs_iBar.end(); index++){ + + for(std::map::iterator index = c_en511_vs_bar.begin(); index != c_en511_vs_bar.end(); index++){ index->second-> Print(Form("%s/summaryPlots/energy/c_energy511_vs_bar_Vov%.01f_th%d.png",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); index->second-> Print(Form("%s/summaryPlots/energy/c_energy511_vs_bar_Vov%.01f_th%d.pdf",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); } - - for(std::map::iterator index = c_en511_vs_iBar_L.begin(); index != c_en511_vs_iBar_L.end(); index++){ - index->second-> Print(Form("%s/summaryPlots/energy/c_energy511_vs_barL_Vov%.01f_th%d.png",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); - index->second-> Print(Form("%s/summaryPlots/energy/c_energy511_vs_barL_Vov%.01f_th%d.pdf",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); - } - - for(std::map::iterator index = c_en511_vs_iBar_R.begin(); index != c_en511_vs_iBar_R.end(); index++){ - index->second-> Print(Form("%s/summaryPlots/energy/c_energy511_vs_barR_Vov%.01f_th%d.png",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); - index->second-> Print(Form("%s/summaryPlots/energy/c_energy511_vs_barR_Vov%.01f_th%d.pdf",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); - } - - - - for(std::map::iterator index = c_en1275_vs_iBar.begin(); index != c_en1275_vs_iBar.end(); index++){ + + for(std::map::iterator index = c_en1275_vs_bar.begin(); index != c_en1275_vs_bar.end(); index++){ index->second-> Print(Form("%s/summaryPlots/energy/c_energy1275_vs_bar_Vov%.01f_th%d.png",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); index->second-> Print(Form("%s/summaryPlots/energy/c_energy1275_vs_bar_Vov%.01f_th%d.pdf",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); } - - for(std::map::iterator index = c_en1275_vs_iBar_L.begin(); index != c_en1275_vs_iBar_L.end(); index++){ - index->second-> Print(Form("%s/summaryPlots/energy/c_energy1275_vs_barL_Vov%.01f_th%d.png",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); - index->second-> Print(Form("%s/summaryPlots/energy/c_energy1275_vs_barL_Vov%.01f_th%d.pdf",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); - } - - for(std::map::iterator index = c_en1275_vs_iBar_R.begin(); index != c_en1275_vs_iBar_R.end(); index++){ - index->second-> Print(Form("%s/summaryPlots/energy/c_energy1275_vs_barR_Vov%.01f_th%d.png",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); - index->second-> Print(Form("%s/summaryPlots/energy/c_energy1275_vs_barR_Vov%.01f_th%d.pdf",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); - } - - - for(std::map::iterator index = c_enRatio_vs_iBar.begin(); index != c_enRatio_vs_iBar.end(); index++){ + + for(std::map::iterator index = c_enRatio_vs_bar.begin(); index != c_enRatio_vs_bar.end(); index++){ index->second-> Print(Form("%s/summaryPlots/energy/c_energyRatio_vs_bar_Vov%.01f_th%d.png",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); index->second-> Print(Form("%s/summaryPlots/energy/c_energyRatio_vs_bar_Vov%.01f_th%d.pdf",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); } - - for(std::map::iterator index = c_enRatio_vs_iBar_L.begin(); index != c_enRatio_vs_iBar_L.end(); index++){ - index->second-> Print(Form("%s/summaryPlots/energy/c_energyRatiovs_barL_Vov%.01f_th%d.png",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); - index->second-> Print(Form("%s/summaryPlots/energy/c_energyRatio_vs_barL_Vov%.01f_th%d.pdf",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); - } - - for(std::map::iterator index = c_enRatio_vs_iBar_R.begin(); index != c_enRatio_vs_iBar_R.end(); index++){ - index->second-> Print(Form("%s/summaryPlots/energy/c_energyRatio_vs_barR_Vov%.01f_th%d.png",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); - index->second-> Print(Form("%s/summaryPlots/energy/c_energyRatio_vs_barR_Vov%.01f_th%d.pdf",plotDir.c_str(),float((int(index->first/10000))/100),int((index->first - (float((int(index->first/10000))/100))*1000000)/100))); - } - - - -//Time Resolution - + + + + + // Time Resolution std::map map_timeRes; std::map map_timeRes511; std::map map_timeRes1275; @@ -1694,8 +1594,8 @@ int main(int argc, char** argv) std::map g_timeRes_vs_th; std::map g_timeRes_vs_Vov; - std::map g_timeRes511_vs_iBar; - std::map g_timeRes1275_vs_iBar; + std::map g_timeRes511_vs_bar; + std::map g_timeRes1275_vs_bar; std::map g_timeRes_vs_enBin; for (auto mapIt : trees2){ @@ -1762,14 +1662,14 @@ int main(int argc, char** argv) } if (enBin ==3){ - if( g_timeRes511_vs_iBar[Vov_th_ID] == NULL ){ - g_timeRes511_vs_iBar[Vov_th_ID] = new TGraphErrors(); + if( g_timeRes511_vs_bar[Vov_th_ID] == NULL ){ + g_timeRes511_vs_bar[Vov_th_ID] = new TGraphErrors(); } } if (enBin ==8){ - if( g_timeRes1275_vs_iBar[Vov_th_ID] == NULL ){ - g_timeRes1275_vs_iBar[Vov_th_ID] = new TGraphErrors(); + if( g_timeRes1275_vs_bar[Vov_th_ID] == NULL ){ + g_timeRes1275_vs_bar[Vov_th_ID] = new TGraphErrors(); } } @@ -1789,15 +1689,15 @@ int main(int argc, char** argv) } if (enBin ==3){ - g_timeRes511_vs_iBar[Vov_th_ID]->SetPoint(g_timeRes511_vs_iBar[Vov_th_ID]->GetN(), iBar, (index->second)/2); - //std::cout<<"N"<GetN()<<"iBar"<second)/2<SetPointError(g_timeRes511_vs_iBar[Vov_th_ID]->GetN()-1, 0., 0.); + g_timeRes511_vs_bar[Vov_th_ID]->SetPoint(g_timeRes511_vs_bar[Vov_th_ID]->GetN(), iBar, (index->second)/2); + //std::cout<<"N"<GetN()<<"iBar"<second)/2<SetPointError(g_timeRes511_vs_bar[Vov_th_ID]->GetN()-1, 0., 0.); } if (enBin ==8){ - g_timeRes1275_vs_iBar[Vov_th_ID]->SetPoint(g_timeRes1275_vs_iBar[Vov_th_ID]->GetN(), iBar, (index->second)/2); - //std::cout<<"N"<GetN()<<"iBar"<second)/2<SetPointError(g_timeRes1275_vs_iBar[Vov_th_ID]->GetN()-1, 0., 0.); + g_timeRes1275_vs_bar[Vov_th_ID]->SetPoint(g_timeRes1275_vs_bar[Vov_th_ID]->GetN(), iBar, (index->second)/2); + //std::cout<<"N"<GetN()<<"iBar"<second)/2<SetPointError(g_timeRes1275_vs_bar[Vov_th_ID]->GetN()-1, 0., 0.); } @@ -1810,7 +1710,7 @@ int main(int argc, char** argv) std::map c_timeRes_vs_th; std::map c_timeRes_vs_Vov; std::map c_timeRes_vs_enBin; - std::map c_timeRes_vs_iBar; + std::map c_timeRes_vs_bar; std::map c_timeRes_vs_enBin_thBest; std::map tBest; double Vov_iBar_enBin_ID; @@ -1852,6 +1752,8 @@ int main(int argc, char** argv) g_tRes -> SetMarkerColor(1+iter[iBar_enBin_ID]); g_tRes -> SetMarkerStyle(20); g_tRes -> Draw("PL,same"); + outFile -> cd(); + g_tRes -> Write(Form("g_timeRes_vs_th_Vov%.01f_bar%02d_enBin%d",Vov,bar,energyBin)); std::string VovLabel = Form("Vov%.01f", Vov); latex = new TLatex(0.55,0.85-0.04*iter[iBar_enBin_ID],VovLabel.c_str()); @@ -1921,6 +1823,9 @@ int main(int argc, char** argv) g_tRes -> SetMarkerColor(1+iter[iBar_enBin_ID]); g_tRes -> SetMarkerStyle(20); g_tRes -> Draw("PL,same"); + outFile -> cd(); + g_tRes -> Write(Form("g_timeRes_vs_Vov_th%d_bar%02d_enBin%d",th,bar,energyBin)); + std::string thLabel = Form("th%d", th); latex = new TLatex(0.55,0.85-0.04*iter[iBar_enBin_ID],thLabel.c_str()); @@ -1990,12 +1895,16 @@ int main(int argc, char** argv) g_tRes_en -> SetMarkerColor(kBlack); g_tRes_en -> SetMarkerStyle(20); g_tRes_en -> Draw("PL"); + outFile -> cd(); + g_tRes_en -> Write(Form("g_timeRes_vs_enBin_bar%02d_Vov%.01f_th%d",iBar,Vov, th)); + + for ( int nPoint = 0; nPoint < g_timeRes_vs_enBin[Vov_th_iBar_ID]->GetN(); nPoint++){ float eBin = g_timeRes_vs_enBin[Vov_th_iBar_ID]->GetPointX(nPoint); - if(vec_th[i] == tBest[double(10000000*(nPoint+1)) + double(vec_Vov[i]*1000000) + double(bar)]){ + if(vec_th[i] == tBest[double(10000000*(nPoint+1)) + double(Vov*1000000) + double(bar)]){ vec_tRes[iBar_Vov_ID].push_back(std::make_pair ( eBin, g_timeRes_vs_enBin[Vov_th_iBar_ID]->GetPointY(nPoint))); vec_tRes_error[iBar_Vov_ID].push_back(std::make_pair ( g_timeRes_vs_enBin[Vov_th_iBar_ID]->GetErrorX(nPoint), g_timeRes_vs_enBin[Vov_th_iBar_ID]->GetErrorY(nPoint))); @@ -2044,6 +1953,9 @@ int main(int argc, char** argv) graph-> SetMarkerColor(kBlack); graph -> SetMarkerStyle(20); graph-> Draw("PL"); + outFile -> cd(); + graph -> Write(Form("g_timeRes_vs_enBin_bar%02d_Vov%.01f_bestTh",iBar,Vov)); + index->second-> Print(Form("%s/summaryPlots/timeResolution/c_timeRes_vs_enBin_bar%02d_Vov%.01f_bestTh.png",plotDir.c_str(),iBar,Vov)); index->second-> Print(Form("%s/summaryPlots/timeResolution/c_timeRes_vs_enBin_bar%02d_Vov%.01f_bestTh.pdf",plotDir.c_str(),iBar,Vov)); } @@ -2054,7 +1966,7 @@ int main(int argc, char** argv) std::map>>vec_tRes_error_511; std::map>>vec_tRes_error_1275; std::map>>vec_tRes_1275; - for(std::map::iterator index = g_timeRes511_vs_iBar.begin(); index != g_timeRes511_vs_iBar.end(); index++){ + for(std::map::iterator index = g_timeRes511_vs_bar.begin(); index != g_timeRes511_vs_bar.end(); index++){ Vov = float(int(index->first/10000)/100.); th = int((index->first - Vov*1000000)/100); @@ -2062,9 +1974,9 @@ int main(int argc, char** argv) Vov_th_ID = Vov*1000000 + th*100; - if(c_timeRes_vs_iBar[Vov]==NULL){ - c_timeRes_vs_iBar[Vov] = new TCanvas(Form("c_timeRes_vs_bar_Vov%.01f",Vov),Form("c_timeRes_vs_bar_Vov%.01f",Vov)); - c_timeRes_vs_iBar[Vov]->cd(); + if(c_timeRes_vs_bar[Vov]==NULL){ + c_timeRes_vs_bar[Vov] = new TCanvas(Form("c_timeRes_vs_bar_Vov%.01f",Vov),Form("c_timeRes_vs_bar_Vov%.01f",Vov)); + c_timeRes_vs_bar[Vov]->cd(); TH1F* hPad = (TH1F*)( gPad->DrawFrame(0.,0.,17.,500.) ); hPad -> SetTitle(";ID bar;#sigma_{t_{Diff}}/2 [ps]"); hPad -> Draw(); @@ -2072,8 +1984,8 @@ int main(int argc, char** argv) } - TGraph* g_timeRes511 = g_timeRes511_vs_iBar[Vov_th_ID]; - TGraph* g_timeRes1275 = g_timeRes1275_vs_iBar[Vov_th_ID]; + TGraph* g_timeRes511 = g_timeRes511_vs_bar[Vov_th_ID]; + TGraph* g_timeRes1275 = g_timeRes1275_vs_bar[Vov_th_ID]; for ( int nPoint = 0; nPoint < g_timeRes511->GetN(); nPoint++){ if (th == tBest[ double(10000000*3) + double(Vov*1000000) + double(int(g_timeRes511->GetPointX(nPoint)))]){ @@ -2093,7 +2005,7 @@ int main(int argc, char** argv) } } - for(std::map::iterator index = c_timeRes_vs_iBar.begin(); index != c_timeRes_vs_iBar.end(); index++){ + for(std::map::iterator index = c_timeRes_vs_bar.begin(); index != c_timeRes_vs_bar.end(); index++){ index->second -> cd(); //Vov = index->first; @@ -2119,6 +2031,9 @@ int main(int argc, char** argv) graph1-> SetMarkerColor(kRed); graph1 -> SetMarkerStyle(20); graph1-> Draw("PL,same"); + outFile -> cd(); + graph1 -> Write(Form("g_timeRes511_vs_bar_Vov%.01f_bestTh",Vov)); + int up2 = vec_tRes_1275[index->first].size(); for (int k = 0; k < up2; k++){ @@ -2139,6 +2054,9 @@ int main(int argc, char** argv) graph2-> SetMarkerColor(kBlue); graph2 -> SetMarkerStyle(20); graph2-> Draw("PL, same"); + outFile -> cd(); + graph2 -> Write(Form("g_timeRes1275_vs_bar_Vov%.01f_bestTh",Vov)); + std::string Label1 = "511 keV Peak"; TLatex* latex1 = new TLatex(0.55,0.85,Label1.c_str()); @@ -2160,6 +2078,16 @@ int main(int argc, char** argv) } + +int bytes = outFile -> Write(); + std::cout << "============================================" << std::endl; + std::cout << "nr of B written: " << int(bytes) << std::endl; + std::cout << "nr of KB written: " << int(bytes/1024.) << std::endl; + std::cout << "nr of MB written: " << int(bytes/1024./1024.) << std::endl; + std::cout << "============================================" << std::endl; + +outFile -> Close(); + } diff --git a/scripts/submit_moduleCharacterization.py b/scripts/submit_moduleCharacterization.py index 509d37d1..93044a62 100755 --- a/scripts/submit_moduleCharacterization.py +++ b/scripts/submit_moduleCharacterization.py @@ -9,91 +9,122 @@ parser = argparse.ArgumentParser(description='This script splits moduleCharacterization_step* tasks in multiple parallel jobs') parser.add_argument("-l", "--label", required=True, type=str, help="job label") +parser.add_argument("-i", "--inputFolder", required=True, type=str, help="input file folder") +parser.add_argument("-o", "--outputFolder", required=True, type=str, help="output file folder") parser.add_argument("-b", "--baseFolder", required=True, type=str, help="base folder") parser.add_argument("-e", "--exeName", required=True, type=str, help="absolute path of executable") parser.add_argument("-r", "--runs", required=True, type=str, help="comma-separated list of runs to be processed") parser.add_argument("-c", "--configFile", required=True, type=str, help="config file") -parser.add_argument("-s", "--submit", help="submit jobs", action='store_true') -parser.add_argument("-v", "--verbose", help="increase output verbosity", action='store_true') +parser.add_argument("-s", "--submit", help="submit jobs", action='store_true') +parser.add_argument("-v", "--verbose", help="increase output verbosity", action='store_true') +parser.add_argument("-p", "--numParal", required=True, type=int, help='number of parallel jobs, 0 for all') args = parser.parse_args() +preruns = [] runs = [] comma_list = args.runs.split(',') for item in comma_list: - hyphen_list = item.split('-') - if len(hyphen_list) > 1: - for i in range(int(hyphen_list[0]), int(hyphen_list[1])+1): - runs.append(i) - else: - runs.append(int(hyphen_list[0])) + hyphen_list = item.split('-') + if len(hyphen_list) > 1: + for i in range(int(hyphen_list[0]), int(hyphen_list[1])+1): + preruns.append(i) + else: + preruns.append(int(hyphen_list[0])) + +if args.numParal == 0: + listRuns= [] + for i in range (len(preruns)): + listRuns.append(int(preruns[i])) + runs.append(listRuns); + +if args.numParal !=0: + step = (int)(len(preruns)/args.numParal) + + for i in range (step): + listRuns= [] + for j in range ( args.numParal * i , args.numParal * (i+1)): + listRuns.append(int(preruns[j])) + runs.append(listRuns) + del listRuns + + if len(preruns) % args.numParal != 0: + listRuns= [] + for j in range ( step * args.numParal, len(preruns)): + listRuns.append(int(preruns[j])) + runs.append(listRuns) + + print runs - print print 'START' print -currDir = os.getcwd() -print -try: - subprocess.check_output(['mkdir','jobs']) -except subprocess.CalledProcessError as e: - print e.output -try: - subprocess.check_output(['mkdir','jobs/'+args.label]) -except subprocess.CalledProcessError as e: - print e.output - - -parallelCommand = "parallel --results ./jobs/"+args.label+"/ "+args.baseFolder+"/"+args.exeName+" ::: " - - -##### creates directory and file list for job ####### -jobDir = currDir+'/jobs/'+args.label+'/' -os.system('mkdir '+jobDir) -os.chdir(jobDir) - -for run in runs: - ##### creates config file ####### - with open(args.baseFolder+'/'+args.configFile) as fi: - contents = fi.read() - configFileName = jobDir+"/config_run"+str(run)+".cfg" - with open(configFileName, "w") as fo: - fo.write(contents) - command = 'sed -i \"s%^runs .*$%runs '+str(run)+'%\" '+configFileName - os.system(command) - command = 'sed -i \"s%^step1FileName .*$%step1FileName '+args.baseFolder+'/plots/moduleCharacterization_step1_run'+str(run)+'.root%\" '+configFileName - os.system(command) - command = 'sed -i \"s%^outFileNameStep1 .*$%outFileNameStep1 '+args.baseFolder+'/plots/moduleCharacterization_step1_run'+str(run)+'.root%\" '+configFileName - os.system(command) - command = 'sed -i \"s%^outFileNameStep2 .*$%outFileNameStep2 '+args.baseFolder+'/plots/moduleCharacterization_step2_run'+str(run)+'.root%\" '+configFileName - os.system(command) +for i in range (len(runs)): + currDir = os.getcwd() + + print + + try: + subprocess.check_output(['mkdir','jobs']) + except subprocess.CalledProcessError as e: + print e.output + try: + subprocess.check_output(['mkdir','jobs/'+args.label]) + except subprocess.CalledProcessError as e: + print e.output + + + parallelCommand = "parallel --results ./jobs/"+args.label+"/ "+args.baseFolder+"/"+args.exeName+" ::: " + + + ##### creates directory and file list for job ####### + jobDir = args.baseFolder+'/jobs/'+args.label+'/' + os.system('mkdir '+jobDir) + os.chdir(jobDir) + for run in runs[i]: + ##### creates config file ####### + with open(args.baseFolder + '/'+args.configFile) as fi: + contents = fi.read() + configFileName = jobDir+"/config_run"+str(run)+".cfg" + with open(configFileName, "w") as fo: + fo.write(contents) + command = 'sed -i \"s%^runs .*$%runs '+str(run)+'%\" '+configFileName + os.system(command) + command = 'sed -i \"s%^step1FileName .*$%step1FileName '+args.inputFolder+'/plots/moduleCharacterization_step1_run'+str(run)+'.root%\" '+configFileName + os.system(command) + command = 'sed -i \"s%^outFileNameStep1 .*$%outFileNameStep1 '+args.outputFolder+'/plots/moduleCharacterization_step1_run'+str(run)+'.root%\" '+configFileName + os.system(command) + command = 'sed -i \"s%^outFileNameStep2 .*$%outFileNameStep2 '+args.outputFolder+'/plots/moduleCharacterization_step2_run'+str(run)+'.root%\" '+configFileName + os.system(command) + #command = 'sed -i \"s%^plotDir .*$%plotDir '+'/var/www/html/ModuleCharacterization/ALIO/Co60/FBK_thinQuartz/Data/'+str(run)+ '%\" '+configFileName + #os.system(command) - parallelCommand += configFileName + " " - - -##### creates job file ####### -with open(jobDir+'/jobs.sh', 'w') as fout: - fout.write("#!/bin/sh\n") - fout.write("echo\n") - fout.write("echo 'START---------------'\n") - fout.write("echo 'current dir: ' ${PWD}\n") - fout.write("cd "+str(args.baseFolder)+"\n") - fout.write("echo 'current dir: ' ${PWD}\n") - fout.write("source scripts/setup.sh\n") - fout.write(parallelCommand+"\n") - fout.write("echo 'STOP---------------'\n") - fout.write("echo\n") - fout.write("echo\n") - os.system("chmod 755 job.sh") + parallelCommand += configFileName + " " + + + ##### creates job file ####### + with open(jobDir+'/jobs.sh', 'w') as fout: + fout.write("#!/bin/sh\n") + fout.write("echo\n") + fout.write("echo 'START---------------'\n") + fout.write("echo 'current dir: ' ${PWD}\n") + fout.write("cd "+str(args.baseFolder)+"\n") + fout.write("echo 'current dir: ' ${PWD}\n") + fout.write("source scripts/setup.sh\n") + fout.write(parallelCommand+"\n") + fout.write("echo 'STOP---------------'\n") + fout.write("echo\n") + fout.write("echo\n") + os.system("chmod 755 job.sh") -if args.submit: - os.system("source "+jobDir+"/jobs.sh") - + if args.submit: + os.system("source "+jobDir+"/jobs.sh") + print 'end job' print print 'END' print diff --git a/src/Co60SpectrumAnalyzer.cc b/src/Co60SpectrumAnalyzer.cc index 0284ba28..1bfcef61 100644 --- a/src/Co60SpectrumAnalyzer.cc +++ b/src/Co60SpectrumAnalyzer.cc @@ -6,8 +6,7 @@ std::map > Co60SpectrumAnalyzer(TH1F* histo, std::vector* ranges) { gStyle -> SetOptFit(0); - - + histo->GetXaxis()->SetRangeUser(0,50); std::map > res; @@ -17,12 +16,12 @@ std::map > Co60SpectrumAnalyzer(TH1F* histo, while ( histo -> GetBinContent(startBin) < 100 && histo -> GetBinContent(startBin-1) < 100 ){ startBin ++; } - - endBin = startBin; - while ( histo -> GetBinContent(endBin) > 2 && histo -> GetBinContent(endBin-1) > 2 && histo -> GetBinContent(endBin+1) > 2){ + + endBin = histo->GetNbinsX() /2.; + while ( histo -> GetBinContent(endBin) < 50 && histo -> GetBinContent(endBin-1) < 50 && histo -> GetBinContent(endBin+1) < 50 && endBin < histo->GetNbinsX() -50){ endBin++; } - + //Derivative calculation std::vector derivative; derivative.push_back(0); @@ -30,19 +29,17 @@ std::map > Co60SpectrumAnalyzer(TH1F* histo, if ( bin < startBin){ derivative.push_back(0); - } - if( bin >= startBin){ - - double difference = histo->GetBinContent(bin+1) - histo->GetBinContent(bin); - double distance = histo->GetBinCenter(bin+1) - histo->GetBinCenter(bin); - derivative.push_back(difference/distance); - } + } + if( bin >= startBin){ + double difference = histo->GetBinContent(bin+1) - histo->GetBinContent(bin); + double distance = histo->GetBinCenter(bin+1) - histo->GetBinCenter(bin); + derivative.push_back(difference/distance); + } } - + //Looking for the right side of the 1332 Peak std::vector rangesDerivative; double drop = startBin+1; - int counter = 0; double sum = 0; for ( int j = startBin; j < endBin; j++){ @@ -53,152 +50,155 @@ std::map > Co60SpectrumAnalyzer(TH1F* histo, rangesDerivative.push_back(double(sum)/5.); counter = 0; sum = 0; + } } - + int check = 0; - - while ( check != 99){ - double min = *std::min_element( rangesDerivative.begin(), rangesDerivative.end()); + std::vector erasedIndex; + erasedIndex.push_back(0); + erasedIndex.push_back(rangesDerivative.size()); int minIndex = 5; - do{minIndex ++;} - while(rangesDerivative[minIndex] != min); - - check = 99; - - for ( int i = minIndex +1; i < minIndex + 5; i++){ - if( rangesDerivative[i] > 0){ + double min; + while ( check != 99){ + min = 10000; + for( int j = 0; j < erasedIndex.size()-1; j++){ + for ( int i = erasedIndex[j] +1; i < erasedIndex[j+1] ; i++){ + if(rangesDerivative[i] < min){ + min = rangesDerivative[i]; + minIndex = i; + } + } + } + check = 99; + int risalite = 0; + int esclusioni = 0; + for ( int i = minIndex -1; i < minIndex + 4; i++){ + if( rangesDerivative[i] > 0){ + risalite ++; + if( rangesDerivative[i] > 160){ + esclusioni ++; + } + } + } + if( risalite > 1 || esclusioni > 0){ check = 0; } + if (check == 0){ + erasedIndex.push_back(minIndex); + } + sort(erasedIndex.begin(), erasedIndex.end()); + drop = histo -> GetBinCenter(startBin + (minIndex+5) *5); } - if (check == 0){ - rangesDerivative.erase(rangesDerivative.begin() + minIndex); - } - drop = histo -> GetBinCenter(startBin + (minIndex+5) *5); + std::vector totalPeaks; + std::vector realPeaks; - } + TF1* fitFunc; - //Searching 1173 peak - int nPeaks = 4; - TSpectrum* spectrum = new TSpectrum(nPeaks); - - histo -> GetXaxis() -> SetRangeUser(0., drop); - - int nFound = spectrum -> Search(histo, 0.01, "", 0.001); - double* peaks = spectrum -> GetPositionX(); + if (drop > histo->GetBinCenter((endBin-startBin)/3.) ){ + //Searching 1173 and 1332 peaks + int nPeaks = 4; + int nFound = 0; + TSpectrum* spectrum = new TSpectrum(nPeaks); + float start = histo->GetBinCenter(startBin); + float stop = histo->GetBinCenter(endBin); + histo -> GetXaxis() -> SetRangeUser(start, stop); + + nFound = spectrum -> Search(histo, 1,"" , 0.001); + double* peaks = spectrum -> GetPositionX(); + int index = -1; + if (nFound > 1){ + for( int i = 0; i < nFound; i++){ + totalPeaks.push_back(peaks[i]); + } + std::sort(totalPeaks.begin(),totalPeaks.end()); + for( int i=0; i< totalPeaks.size(); i++){ + if( totalPeaks[i] < drop){ + index ++; + } + } + int a; + int b; + if( histo -> GetBinContent( histo -> FindBin(totalPeaks[index])) > histo -> GetBinContent( histo -> FindBin(totalPeaks[index-1]))){ + a = index; + b = index -1; + } + if( histo -> GetBinContent( histo -> FindBin(totalPeaks[index])) < histo -> GetBinContent( histo -> FindBin(totalPeaks[index-1]))){ + b = index; + a = index -1; + } - std::vector totalPeaks; - std::vector totalPeaks2; - std::vector realPeaks; - - - if(nFound >1){ - for( int i = 0; i < nFound; i++){ - totalPeaks.push_back(peaks[i]); - } - - sort(totalPeaks.begin(), totalPeaks.end()); - - for( int j = 0; j < totalPeaks.size(); j++){ - } - } - - if(totalPeaks.size() == 4){ - realPeaks.push_back(totalPeaks[totalPeaks.size()-2]); - realPeaks.push_back(totalPeaks[totalPeaks.size()-1]); - } - if(totalPeaks.size() <4){ - double d1 = drop - histo->GetBinWidth(10)*20 -totalPeaks[totalPeaks.size()-1]; - TF1* fitFunc_gaus = new TF1("fitFunc_gaus","gaus",0.95*totalPeaks[2],1.05*totalPeaks[2]); - histo -> Fit(fitFunc_gaus,"QR"); - double sigma = fitFunc_gaus -> GetParameter(2); - if(d1 > sigma*1.5){ - realPeaks.push_back(totalPeaks[totalPeaks.size()-1]); - } - if(d1 < sigma*1.5){ - realPeaks.push_back(totalPeaks[totalPeaks.size()-2]); - realPeaks.push_back(totalPeaks[totalPeaks.size()-1]); + if( histo -> GetBinContent( histo -> FindBin(totalPeaks[a])) / histo -> GetBinContent( histo -> FindBin(totalPeaks[b])) < 10){ + realPeaks.push_back(totalPeaks[index-1]); + realPeaks.push_back(totalPeaks[index]); + } + if( histo -> GetBinContent( histo -> FindBin(totalPeaks[a])) / histo -> GetBinContent( histo -> FindBin(totalPeaks[b])) > 10){ + realPeaks.push_back(totalPeaks[index-2]); + realPeaks.push_back(totalPeaks[index-1]); + } + + } } - } - - - - //Fitting 1173 keV peak - if (realPeaks.size() == 1 || realPeaks.size() == 2 ){ - TF1* fitFunc_1173 = new TF1("fitFunc_1173","gaus",0.95*realPeaks[0],1.05*realPeaks[0]); - fitFunc_1173 -> SetLineColor(kBlack); - fitFunc_1173 -> SetLineWidth(2); - histo -> Fit(fitFunc_1173,"QRS+"); - - res["1.173 MeV"] = std::make_pair(fitFunc_1173->GetParameter(1),fitFunc_1173->GetParameter(2)); - //histo -> GetXaxis() -> SetRangeUser(0.,3.*fitFunc_1173->GetParameter(1)); - } - - //Searching 1332 peaks if previously failed - if(realPeaks.size() == 1){ - histo -> GetXaxis() -> SetRangeUser ( res["1.173 MeV"].first +res["1.173 MeV"].second, 40.); - nPeaks = 2; - TSpectrum* spectrum2 = new TSpectrum(nPeaks); - nFound = spectrum2 -> Search(histo, 0.01, "", 0.001); - double* peaks2 = spectrum2 -> GetPositionX(); - - if(nFound >0){ - for( int i = 0; i < nFound; i++){ - totalPeaks2.push_back(peaks2[i]); - } - - sort(totalPeaks2.begin(), totalPeaks2.end()); - - realPeaks.push_back(totalPeaks2[0]); - } - } - - histo -> GetXaxis() -> SetRangeUser(0., 40.); - + histo -> GetXaxis() -> SetRangeUser(0, 50); + if (realPeaks.size() == 1 || realPeaks.size() == 2 ){ + fitFunc = new TF1( "fitFunc", "gaus", realPeaks[0] -0.5, realPeaks[0] +0.5); + histo-> Fit(fitFunc, "NQR"); + TF1* fitFunc_1173 = new TF1("fitFunc_1173","gaus",fitFunc->GetParameter(1)-fitFunc->GetParameter(2)*0.75, fitFunc->GetParameter(1)+fitFunc->GetParameter(2)*0.5); + fitFunc_1173 -> SetParameter(0, fitFunc->GetParameter(0)); + fitFunc_1173 -> SetParameter(1, fitFunc->GetParameter(1)); + fitFunc_1173 -> SetParameter(2, fitFunc->GetParameter(2)); + fitFunc_1173 -> SetLineColor(kBlack); + fitFunc_1173 -> SetLineWidth(2); + histo -> Fit(fitFunc_1173,"QRS+"); + res["1.173 MeV"] = std::make_pair(fitFunc_1173->GetParameter(1),fitFunc_1173->GetParameter(2)); + //histo -> GetXaxis() -> SetRangeUser(0.,3.*fitFunc_1173->GetParameter(1)); + } //Fitting 1332 keV peak if (realPeaks.size() == 2){ - - TF1* fitFunc_1332 = new TF1("fitFunc_1332","gaus",0.95*realPeaks[1],1.05*realPeaks[1]); - fitFunc_1332 -> SetLineColor(kBlack); - fitFunc_1332 -> SetLineWidth(2); - histo -> Fit(fitFunc_1332,"QRS+"); + TF1* fitFunc2 = new TF1( "fitFunc2", "gaus", realPeaks[1] -0.5, realPeaks[1] +0.5); + histo-> Fit(fitFunc2, "NQR"); + TF1* fitFunc_1332 = new TF1("fitFunc_1332","gaus",fitFunc2->GetParameter(1)-fitFunc2->GetParameter(2)*0.75, fitFunc2->GetParameter(1)+fitFunc2->GetParameter(2)*0.75); + fitFunc_1332 -> SetParameter(0, fitFunc2->GetParameter(0)); + fitFunc_1332 -> SetParameter(1, fitFunc2->GetParameter(1)); + fitFunc_1332 -> SetParameter(2, fitFunc2->GetParameter(2)); + fitFunc_1332 -> SetLineColor(kBlack); + fitFunc_1332 -> SetLineWidth(2); + histo -> Fit(fitFunc_1332,"QRS+"); - res["1.332 MeV"] = std::make_pair(fitFunc_1332->GetParameter(1),fitFunc_1332->GetParameter(2)); + res["1.332 MeV"] = std::make_pair(fitFunc_1332->GetParameter(1),fitFunc_1332->GetParameter(2)); } - - //Not Co60 spectrum controll - if (realPeaks.size() == 0 || realPeaks.size() > 2 ){ - std::cout << "Errore" << std::endl; - res["1.173 MeV"] = std::make_pair(-9999,0); - } + //Locating and drawing ranges of interest if ( realPeaks.size() > 0 && realPeaks.size() < 3 ){ - ranges->push_back(res["1.173 MeV"].first-8.*res["1.173 MeV"].second); - ranges->push_back(res["1.173 MeV"].first-5.*res["1.173 MeV"].second); - ranges->push_back(res["1.173 MeV"].first-2.*res["1.173 MeV"].second); - ranges->push_back(res["1.173 MeV"].first+2.*res["1.173 MeV"].second); - ranges->push_back(res["1.173 MeV"].first+5.*res["1.173 MeV"].second); + ranges->push_back(res["1.173 MeV"].first-12.*res["1.173 MeV"].second); + ranges->push_back(res["1.173 MeV"].first-10.5*res["1.173 MeV"].second); + ranges->push_back(res["1.173 MeV"].first-9.*res["1.173 MeV"].second); + ranges->push_back(res["1.173 MeV"].first-7.5*res["1.173 MeV"].second); + ranges->push_back(res["1.173 MeV"].first-6.*res["1.173 MeV"].second); + ranges->push_back(res["1.173 MeV"].first-4.5*res["1.173 MeV"].second); + ranges->push_back(res["1.173 MeV"].first-3.*res["1.173 MeV"].second); + ranges->push_back(res["1.173 MeV"].first-1.5*res["1.173 MeV"].second); } if ( realPeaks.size() == 2 ){ - ranges->push_back(0.5*(res["1.173 MeV"].first+5.*res["1.173 MeV"].second+res["1.332 MeV"].first-5.*res["1.332 MeV"].second)); - ranges->push_back(res["1.332 MeV"].first-5.*res["1.332 MeV"].second); - ranges->push_back(res["1.332 MeV"].first-2.*res["1.332 MeV"].second); - ranges->push_back(res["1.332 MeV"].first+2.*res["1.332 MeV"].second); - } - - - for(auto range: (*ranges)){ - TLine* line = new TLine(range,3.,range,histo->GetBinContent(histo->FindBin(range))); - line -> SetLineWidth(1); - line -> SetLineStyle(7); - line -> Draw("same"); - - - + ranges->push_back(0.5*(res["1.173 MeV"].first+4.*res["1.173 MeV"].second+res["1.332 MeV"].first-4.*res["1.332 MeV"].second)); + ranges->push_back(res["1.332 MeV"].first+1.5*res["1.332 MeV"].second); + ranges->push_back(res["1.332 MeV"].first+3.*res["1.332 MeV"].second); } - + //Not Co60 spectrum controll + if (realPeaks.size() == 0 || realPeaks.size() > 2 ){ + std::cout << "Errore" << std::endl; + res["1.173 MeV"] = std::make_pair(-9999,0); + } + if( realPeaks.size() == 1 || realPeaks.size() == 2){ + for(auto range: (*ranges)){ + TLine* line = new TLine(range,3.,range,histo->GetBinContent(histo->FindBin(range))); + line -> SetLineWidth(1); + line -> SetLineStyle(7); + line -> Draw("same"); + } + } return res; } diff --git a/src/Na22SpectrumAnalyzer.cc b/src/Na22SpectrumAnalyzer.cc index dccacbaf..6b34b8af 100644 --- a/src/Na22SpectrumAnalyzer.cc +++ b/src/Na22SpectrumAnalyzer.cc @@ -7,192 +7,263 @@ std::map > Na22SpectrumAnalyzer(TH1F* histo, { gStyle -> SetOptFit(0); - std::map > res; + + // use TSpectrum to localize peaks + histo -> GetXaxis() -> SetRangeUser(1.,40.); + int nPeaks = 5; TSpectrum* spectrum = new TSpectrum(nPeaks); - - histo -> GetXaxis() -> SetRangeUser(0.,40.); - int nFound = spectrum -> Search(histo, 0.5, "", 0.001); double* peaks = spectrum -> GetPositionX(); - - //totalPeaks -> all found Peaks, nBins -> bin number of found Peaks, realPeaks -> 511 e 1275 + + + // totalPeaks -> all found Peaks, nBins -> bin number of found Peaks, realPeaks -> 511 and 1275 std::vector totalPeaks; std::vector realPeaks; std::vector nBins; - - if(nFound >2){ - - for( int i = 0; i < nFound; i++){ - nBins.push_back(histo -> FindBin(peaks[i])); - totalPeaks.push_back(peaks[i]); - } - sort(nBins.begin(), nBins.end()); - sort(totalPeaks.begin(), totalPeaks.end()); - - //Derivative calculation - std::vector derivative; - derivative.push_back(0); - - int unemptyBin = 2; - while ( histo -> GetBinContent(unemptyBin) < 100 && histo -> GetBinContent(unemptyBin-1) < 100 ){ - unemptyBin ++; - } - for( int bin = 1; bin < nBins[nBins.size()-1] +20; bin++){ - if ( bin < unemptyBin){ - derivative.push_back(0); - } - if( bin >= unemptyBin){ - double difference = histo->GetBinContent(bin+1) - histo->GetBinContent(bin); - double distance = histo->GetBinCenter(bin+1) - histo->GetBinCenter(bin); - derivative.push_back(difference/distance); - } - } - - - - //Finding first and second maximum in totalPeaks - //Fitting them with gaussian and getting their sigma double maxValue = 0; - double secondMaxValue = 0; - int maxIndex = 0; - int secondMaxIndex = 0; - double maxSigma =0; - double secondMaxSigma =0; + double secondMaxValue = 0; + int maxIndex = 0; + int secondMaxIndex = 0; + double maxSigma = 0; + double secondMaxSigma = 0; + int endBin = 0; + int nFoundAdd = 0; + + if(nFound > 1){ + for(int i = 0; i < nFound; ++i) + { + nBins.push_back( histo->FindBin(peaks[i]) ); + totalPeaks.push_back( peaks[i] ); + } + sort(nBins.begin(),nBins.end()); + sort(totalPeaks.begin(),totalPeaks.end()); - //First max - for (int i = 0; i < nFound; i++){ - if ( histo -> GetBinContent(nBins[i]) > maxValue){ - maxValue = histo -> GetBinContent(nBins[i]); - maxIndex = i; + endBin = nBins[0]; + while ( histo -> GetBinContent(endBin) > 10 && histo -> GetBinContent(endBin-1) > 10 && histo -> GetBinContent(endBin+1) > 10 ){ + endBin++; } - } - TF1* fitFunc_gaus = new TF1("fitFunc_gaus","gaus",0.95*double(totalPeaks[maxIndex]),1.05*double(totalPeaks[maxIndex])); - histo -> Fit(fitFunc_gaus,"QR"); - maxSigma = fitFunc_gaus->GetParameter(2); - - //Second Max - for (int i = 0; i < nFound; i++){ - if(i != maxIndex){ - if ( histo -> GetBinContent(nBins[i]) > secondMaxValue){ - secondMaxValue = histo -> GetBinContent(nBins[i]); - secondMaxIndex = i; + //looking for one other peak after the found ones + if(totalPeaks[totalPeaks.size()-1]+2 < histo->GetBinCenter(endBin)){ + histo -> GetXaxis() -> SetRangeUser(totalPeaks[totalPeaks.size()-1]*1.1 ,40.); + int nPeaksAdd = 1; + TSpectrum* spectrumAdd = new TSpectrum(nPeaksAdd); + nFoundAdd = spectrumAdd -> Search(histo, 0.5, "", 0.001); + double* peaksAdd = spectrumAdd -> GetPositionX(); + if (nFoundAdd != 0 && peaksAdd[0] != totalPeaks[totalPeaks.size()-1] && histo->FindBin(peaksAdd[0]) < endBin ){ + totalPeaks.push_back(peaksAdd[0]); + nBins.push_back(histo->FindBin(peaksAdd[0])); } } - } - TF1* fitFunc_gaus2 = new TF1("fitFunc_gaus2","gaus",0.95*double(totalPeaks[secondMaxIndex]),1.05*double(totalPeaks[secondMaxIndex])); - histo -> Fit(fitFunc_gaus2,"QR"); - secondMaxSigma = fitFunc_gaus2->GetParameter(2); - - - //Finding compton edge, looking for the range (width 10 bins) with smallest derivative not included in [first max - second max] - int sigmaBin = std::min(maxSigma, secondMaxSigma) / (histo->GetBinWidth(10)); - int windowBin = 10; - double minPendence = 1000; - int compton; - - for ( int j = 1; j < nFound-1; j++){ + + histo -> GetXaxis() -> SetRangeUser(0. ,40.); - int nRanges = 0; - if ( double(nBins[j+1] -sigmaBin * 2 - nBins[j])/windowBin - int((nBins[j+1] -sigmaBin * 2 - nBins[j])/windowBin) < 0.5){ - nRanges = int((nBins[j+1] -sigmaBin * 2 - nBins[j])/windowBin); + + + + if( nFound + nFoundAdd > 2 ) + { + + + //Derivative calculation + std::vector derivative; + derivative.push_back(0); + + int unemptyBin = 2; + while( histo -> GetBinContent(unemptyBin) < 100 && histo -> GetBinContent(unemptyBin-1) < 100 ) + { + ++unemptyBin; + } + for(int bin = 1; bin < nBins[nBins.size()-1]+20; ++bin) + { + if ( bin < unemptyBin) derivative.push_back(0); + else + { + double difference = histo->GetBinContent(bin+1) - histo->GetBinContent(bin); + double distance = histo->GetBinCenter(bin+1) - histo->GetBinCenter(bin); + derivative.push_back(difference/distance); + } + } + + // Finding first and second maximum in totalPeaks + // Fitting them with gaussian and getting their sigma + + //First max + for (int i = 0; i < nFound; i++){ + if ( histo -> GetBinContent(nBins[i]) > maxValue){ + maxValue = histo -> GetBinContent(nBins[i]); + maxIndex = i; + } + } + TF1* fitFunc_gaus = new TF1("fitFunc_gaus","gaus",0.95*double(totalPeaks[maxIndex]),1.05*double(totalPeaks[maxIndex])); + histo -> Fit(fitFunc_gaus,"QR"); + maxSigma = fitFunc_gaus->GetParameter(2); + + //Second Max + if(nFound < 4){ + secondMaxValue = maxValue; + secondMaxIndex = maxIndex; } - if ( double(nBins[j+1] -sigmaBin * 2 - nBins[j])/windowBin - int((nBins[j+1] -sigmaBin * 2. - nBins[j])/windowBin) >= 0.5){ - nRanges = int((nBins[j+1] -sigmaBin * 2 - nBins[j])/windowBin) +1; + if(nFound > 3){ + for (int i = 0; i < nFound; i++){ + if(i != maxIndex){ + if ( histo -> GetBinContent(nBins[i]) > secondMaxValue){ + secondMaxValue = histo -> GetBinContent(nBins[i]); + secondMaxIndex = i; + } + } + } } + TF1* fitFunc_gaus2 = new TF1("fitFunc_gaus2","gaus",0.95*double(totalPeaks[secondMaxIndex]),1.05*double(totalPeaks[secondMaxIndex])); + histo -> Fit(fitFunc_gaus2,"QR"); + secondMaxSigma = fitFunc_gaus2->GetParameter(2); + + //Finding compton edge between 511 e 1275 + int sigmaBin = std::min(maxSigma, secondMaxSigma) / (histo->GetBinWidth(10)); + int windowBin = 5; + double minPendence = 1000; + int compton = totalPeaks.size()-1; + + for ( int j = 0; j < nFound-1; j++) + { + int nRanges = 0; + if ( double(nBins[j+1] -sigmaBin * 2 - nBins[j])/windowBin - int((nBins[j+1] -sigmaBin * 2 - nBins[j])/windowBin) < 0.5){ + nRanges = int((nBins[j+1] -sigmaBin * 2 - nBins[j])/windowBin); + } + if ( double(nBins[j+1] -sigmaBin * 2 - nBins[j])/windowBin - int((nBins[j+1] -sigmaBin * 2. - nBins[j])/windowBin) >= 0.5){ + nRanges = int((nBins[j+1] -sigmaBin * 2 - nBins[j])/windowBin) +1; + } + + for(int r = 0; r < nRanges; r++){ + double pendence = 0; + int counter = 0; + for ( int i = nBins[j]+ sigmaBin + windowBin*r; i < nBins[j] + sigmaBin +windowBin*(r+1); i++){ + pendence += derivative[i]; + counter ++; + } - for(int r = 0; r < nRanges; r++){ - double pendence = 0; - for ( int i = nBins[j]+ sigmaBin + windowBin*r; i < nBins[j] + sigmaBin +windowBin*(r+1); i++){ - pendence += derivative[i]; + if(std::abs(pendence/counter) <= minPendence && j+1 > maxIndex && j+1 > secondMaxIndex && histo->GetBinContent(nBins[j+1]) > 50 ){ //std::abs(pendence/counter)>0.1 && + minPendence = std::abs(pendence/counter); + compton = j+1; + } } - - if(std::abs(pendence/windowBin) <= minPendence && std::abs(pendence/windowBin)>0.1 && j+1 > maxIndex && j+1 > secondMaxIndex){ - minPendence = std::abs(pendence/windowBin); - compton = j+1; + } + int afterCompton = nFound + nFoundAdd - 1 - compton; + // Compton is the i-th found peaks after the compton edge + // Looking for 511 peak valuating the amplitude of peaks before compton + + int firstPeak = 0; + int found1275 = 0; + while ( found1275 == 0){ + firstPeak++; + if ( double(histo->GetBinContent(nBins[compton-firstPeak]))/double(secondMaxValue) > 0.8 || double(histo->GetBinContent(nBins[compton-firstPeak]))/double(maxValue) > 0.2 ){ + realPeaks.push_back(totalPeaks[compton-firstPeak]); + found1275 = 1; } } + + //Fitting 511 keV peak + if( realPeaks.size() >= 1 ) + { + TF1*fitFunc_511 = new TF1("fitFunc_511","gaus",realPeaks[0] *0.90,realPeaks[0] *1.1); + fitFunc_511 -> SetLineColor(kBlack); + fitFunc_511 -> SetLineWidth(2); + histo -> Fit(fitFunc_511,"QRS+"); + res["0.511 MeV"] = std::make_pair(fitFunc_511->GetParameter(1),fitFunc_511->GetParameter(2)); + histo -> GetXaxis() -> SetRangeUser(0.,3.*fitFunc_511->GetParameter(1)); } - + + //Looking for 1275 keV peak in a window [1.8* 511 position, 4* 511 position] + int secondFound = 0; + std::vector secondPeaksIndex; + + for( int i = 0 ; i < totalPeaks.size(); i++){ + if ( totalPeaks[i] > res["0.511 MeV"].first * 1.8 && totalPeaks[i] < res["0.511 MeV"].first * 4){ + secondFound ++; + secondPeaksIndex.push_back(i); + } + } + if( secondFound == 1){ + realPeaks.push_back(totalPeaks[secondPeaksIndex[0]]); + } + + if( secondPeaksIndex.size() == 2){ + double peaksRatio = histo->GetBinContent(nBins[secondPeaksIndex[0]]) / histo->GetBinContent(nBins[secondPeaksIndex[1]]); + if ( peaksRatio > 3 ) realPeaks.push_back(totalPeaks[secondPeaksIndex[0]]); + if ( peaksRatio < 3 && peaksRatio > 0.3333) realPeaks.push_back(totalPeaks[secondPeaksIndex[1]]); + } - - //Compton is the i-th found peaks after the compton edge - //Based on the number of peaks before and after the compton edge, 511 and 1275 Peaks are selected - if( compton > maxIndex && compton > secondMaxIndex){ - int firstPeak = 0; - int afterCompton = nFound -1 -compton; - - do{ firstPeak++;} - while(double(histo->GetBinContent(nBins[compton-firstPeak]))/double(maxValue) < 0.1); - - realPeaks.push_back(totalPeaks[compton-firstPeak]); - if(afterCompton == 1 || afterCompton == 0){ - realPeaks.push_back(totalPeaks[compton+afterCompton]); + if( secondPeaksIndex.size() == 3){ + realPeaks.push_back(totalPeaks[secondPeaksIndex[1]]); } - - if(afterCompton > 1){ - double peaksRatio = histo->GetBinContent(nBins[nFound-2]) / histo->GetBinContent(nBins[nFound-1]); - if ( peaksRatio > 3 ) { - realPeaks.push_back(totalPeaks[nFound-2]); - } - if ( peaksRatio < 3 && peaksRatio > 0.3333){ - realPeaks.push_back(totalPeaks[nFound-1]); - } + if( secondPeaksIndex.size() == 4){ + realPeaks.push_back(totalPeaks[secondPeaksIndex[2]]); } - } - } - - - - //Fitting 511 keV peak - if (realPeaks.size() == 1 || realPeaks.size() == 2 ){ - TF1* fitFunc_511 = new TF1("fitFunc_511","gaus",0.90*realPeaks[0],1.10*realPeaks[0]); - fitFunc_511 -> SetLineColor(kBlack); - fitFunc_511 -> SetLineWidth(2); - histo -> Fit(fitFunc_511,"QRS+"); - res["0.511 MeV"] = std::make_pair(fitFunc_511->GetParameter(1),fitFunc_511->GetParameter(2)); - histo -> GetXaxis() -> SetRangeUser(0.,3.*fitFunc_511->GetParameter(1)); } + if( realPeaks.size() ==1 ){ + histo -> GetXaxis() -> SetRangeUser(totalPeaks[totalPeaks.size()-1]+2 ,40.); + int nPeaks2 = 1; + TSpectrum* spectrum2 = new TSpectrum(nPeaks2); + int nFound2 = spectrum2 -> Search(histo, 0.5, "", 0.001); + double* peaks2 = spectrum2 -> GetPositionX(); + if (nFound2 != 0){ + realPeaks.push_back(peaks2[0]); + } + } + + histo -> GetXaxis() -> SetRangeUser(0. ,40.); + //Fitting 1275 keV peak - if (realPeaks.size() == 2){ - TF1* fitFunc_1275 = new TF1("fitFunc_1275","gaus",0.95*realPeaks[1],1.05*realPeaks[1]); - fitFunc_1275 -> SetLineColor(kBlack); - fitFunc_1275 -> SetLineWidth(2); - histo -> Fit(fitFunc_1275,"QRS+"); - res["1.275 MeV"] = std::make_pair(fitFunc_1275->GetParameter(1),fitFunc_1275->GetParameter(2)); + if( realPeaks.size() == 2 ) + { + TF1* fitFunc_1275 = new TF1("fitFunc_1275","gaus",realPeaks[1] *0.90, realPeaks[1] *1.10); + fitFunc_1275 -> SetParameter(1, realPeaks[1]); + fitFunc_1275 -> SetParameter(2, res["0.511 Mev"].second/2.); + fitFunc_1275 -> SetLineColor(kBlack); + fitFunc_1275 -> SetLineWidth(2); + histo -> Fit(fitFunc_1275,"QRS+"); + res["1.275 MeV"] = std::make_pair(fitFunc_1275->GetParameter(1),fitFunc_1275->GetParameter(2)); } - + //Not Na22 spectrum controll - if (realPeaks.size() == 0 || realPeaks.size() > 2 || nFound < 3){ - std::cout << "Errore" << std::endl; - res["0.511 MeV"] = std::make_pair(-9999,0); + if( realPeaks.size() == 0 || realPeaks.size() > 2 ) + { + std::cout << "Errore" << std::endl; + res["0.511 MeV"] = std::make_pair(-9999,0); } - + //Locating and drawing ranges of interest if ( realPeaks.size() > 0 && realPeaks.size() < 3 ){ ranges->push_back(res["0.511 MeV"].first-8.*res["0.511 MeV"].second); ranges->push_back(res["0.511 MeV"].first-5.*res["0.511 MeV"].second); - ranges->push_back(res["0.511 MeV"].first-2.*res["0.511 MeV"].second); - ranges->push_back(res["0.511 MeV"].first+2.*res["0.511 MeV"].second); - ranges->push_back(res["0.511 MeV"].first+5.*res["0.511 MeV"].second); + ranges->push_back(res["0.511 MeV"].first-2.*res["0.511 MeV"].second); + ranges->push_back(res["0.511 MeV"].first+2.*res["0.511 MeV"].second); + ranges->push_back(res["0.511 MeV"].first+5.*res["0.511 MeV"].second); + } if ( realPeaks.size() == 2 ){ - ranges->push_back(0.5*(res["0.511 MeV"].first+5.*res["0.511 MeV"].second+res["1.275 MeV"].first-5.*res["1.275 MeV"].second)); - ranges->push_back(res["1.275 MeV"].first-5.*res["1.275 MeV"].second); - ranges->push_back(res["1.275 MeV"].first-2.*res["1.275 MeV"].second); - ranges->push_back(res["1.275 MeV"].first+2.*res["1.275 MeV"].second); + ranges->push_back(0.5*(res["0.511 MeV"].first+5.*res["0.511 MeV"].second+res["1.275 MeV"].first-5.*res["1.275 MeV"].second)); + ranges->push_back(res["1.275 MeV"].first-5.*res["1.275 MeV"].second); + ranges->push_back(res["1.275 MeV"].first-2.*res["1.275 MeV"].second); + ranges->push_back(res["1.275 MeV"].first+2.*res["1.275 MeV"].second); } + for(auto range: (*ranges)){ - TLine* line = new TLine(range,3.,range,histo->GetBinContent(histo->FindBin(range))); - line -> SetLineWidth(1); - line -> SetLineStyle(7); - line -> Draw("same"); + + if ( range < 0 && realPeaks.size() == 1) { //error controll + std::cout << "Errore" << std::endl; + res["0.511 MeV"] = std::make_pair(-9999,0); + } + TLine* line = new TLine(range,3.,range,histo->GetBinContent(histo->FindBin(range))); + line -> SetLineWidth(1); + line -> SetLineStyle(7); + line -> Draw("same"); }