diff --git a/UserTools/EventDisplay/EventDisplay.cpp b/UserTools/EventDisplay/EventDisplay.cpp index 6aba745f0..5f757ae88 100644 --- a/UserTools/EventDisplay/EventDisplay.cpp +++ b/UserTools/EventDisplay/EventDisplay.cpp @@ -120,12 +120,9 @@ bool EventDisplay::Initialise(std::string configfile, DataModel &data){ ifstream file_lappd(filename_lappd); while (!file_lappd.eof()){ file_lappd>>temp_lappd; - active_lappds.push_back(temp_lappd); + active_lappds_user.emplace(temp_lappd,1); } file_lappd.close(); - for (int i_lappd=0;i_lappdGetNumDetectorsInSet("Tank"); n_lappds = geom->GetNumDetectorsInSet("LAPPD"); n_mrd_pmts = geom->GetNumDetectorsInSet("MRD"); n_veto_pmts = geom->GetNumDetectorsInSet("Veto"); - m_data->CStore.Get("detectorkey_to_lappdid",detectorkey_to_lappdid); + m_data->CStore.Get("lappdid_to_detectorkey",lappdid_to_detectorkey); m_data->CStore.Get("channelkey_to_mrdpmtid",channelkey_to_mrdpmtid); m_data->CStore.Get("channelkey_to_faccpmtid",channelkey_to_faccpmtid); - std::cout <<"PID: Num Tank PMTs: "< >* Detectors = geom->GetDetectors(); std::cout <<"Detectors size: "<size()<::iterator it = Detectors->at("LAPPD").begin(); + it != Detectors->at("LAPPD").end(); + ++it){ + Detector *alappd = it->second; + std::cout <<"LAPPD detector key: "<GetDetectorID()<GetDetectorID()); + } + max_num_lappds = lappd_detkeys.size(); + std::cout <<"Number of LAPPDs: "<second; unsigned long detkey = it->first; + pmt_detkeys.push_back(detkey); unsigned long chankey = apmt->GetChannels()->begin()->first; Position position_PMT = apmt->GetDetectorPosition(); if (verbose > 2) std::cout <<"detkey: "< 0) std::cout <<"Properties of the detector configuration: "< 0) std::cout <<"Max y of Tank PMTs: "< vec_zero; - Position zero_position(0.,0.,0.); - vec_zero.push_back(zero_position); - hits_LAPPDs.assign(n_lappds,vec_zero); //--------------------------------------------------------------- //----initialize TApplication in case of displayed graphics------ @@ -222,18 +242,17 @@ bool EventDisplay::Initialise(std::string configfile, DataModel &data){ //--------------------------------------------------------------- if (draw_histograms){ - canvas_pmt = new TCanvas("canvas_pmt","Tank PMT histograms",900,600); - canvas_pmt_supplementary = new TCanvas("canvas_pmt_supplementary","PMT combined",900,600); - canvas_lappd = new TCanvas("canvas_lappd","LAPPD histograms",900,600); + canvas_pmt = new TCanvas("canvas_pmt","Tank PMT histograms",900,600); + canvas_pmt_supplementary = new TCanvas("canvas_pmt_supplementary","PMT combined",900,600); + canvas_lappd = new TCanvas("canvas_lappd","LAPPD histograms",900,600); } canvas_ev_display=new TCanvas("canvas_ev_display","Event Display",900,900); for (int i_lappd=0; i_lappd < max_num_lappds; i_lappd++){ - time_LAPPDs[i_lappd] = nullptr; - charge_LAPPDs[i_lappd] = nullptr; + time_LAPPDs[i_lappd] = nullptr; + charge_LAPPDs[i_lappd] = nullptr; } - return true; } @@ -325,24 +344,14 @@ bool EventDisplay::Execute(){ mrddigitpmtsthisevent.clear(); mrddigitchargesthisevent.clear(); - for (int i_pmt=0;i_pmt temp_vec; - time_lappd.assign(n_lappds,temp_vec); hits_LAPPDs.clear(); - std::vector vec_zero; - Position zero_position(0.,0.,0.); - vec_zero.push_back(zero_position); - hits_LAPPDs.assign(n_lappds,vec_zero); - - for (int i_lappd=0;i_lappdDelete(); if (time_PMTs) time_PMTs->Delete(); if (charge_time_PMTs) charge_time_PMTs->Delete(); - for (int i=0;iDelete(); if (time_LAPPDs[i]) time_LAPPDs[i]->Delete(); } @@ -372,10 +381,11 @@ bool EventDisplay::Execute(){ charge_time_PMTs->GetXaxis()->SetTitle("time [ns]"); charge_time_PMTs->GetYaxis()->SetTitle("charge [p.e.]"); - for (int i=0;i 1) std::cout <<"Loop through MCParticles..."<size(); particlei++){ MCParticle aparticle = mcparticles->at(particlei); @@ -402,6 +412,7 @@ bool EventDisplay::Execute(){ if (aparticle.GetParentPdg() !=0 ) continue; if (aparticle.GetFlag() !=0 ) continue; if ((aparticle.GetPdgCode() == 11 || aparticle.GetPdgCode() == 13)){ //primary particle for Cherenkov tracking should be muon or electron + n_flag0++; draw_vertex_temp = draw_vertex; draw_ring_temp = draw_ring; truevtx = aparticle.GetStartVertex(); @@ -423,6 +434,12 @@ bool EventDisplay::Execute(){ } else continue; } + + //in case there are no primary particles, don't plot vertex and ring + if (n_flag0==0) { + draw_vertex_temp=false; + draw_ring_temp=false; + } //--------------------------------------------------------------- //-------------------Iterate over MCHits ------------------------ @@ -436,9 +453,10 @@ bool EventDisplay::Execute(){ unsigned long chankey = apair.first; if (verbose > 3) std::cout <<"chankey: "<ChannelToDetector(chankey); - int detkey = thistube->GetDetectorID(); + unsigned long detkey = thistube->GetDetectorID(); if (verbose > 3) std::cout <<"detkey: "<GetDetectorElement()=="Tank"){ + hitpmt_detkeys.push_back(detkey); std::vector& Hits = apair.second; int hits_pmt = 0; int wcsim_id; @@ -446,39 +464,41 @@ bool EventDisplay::Execute(){ charge[detkey] += ahit.GetCharge(); time[detkey] += ahit.GetTime(); hits_pmt++; - } + } time[detkey]/=hits_pmt; //use mean time of all hits on one PMT - total_hits_pmts++; + total_hits_pmts++; } - } + } //--------------------------------------------------------------- //-------------------Fill time hists ---------------------------- //--------------------------------------------------------------- - for (int i_pmt=0;i_pmtFill(charge[i_pmt]); - if (time[i_pmt]!=0) time_PMTs->Fill(time[i_pmt]); - if (charge[i_pmt]!=0) charge_time_PMTs->Fill(time[i_pmt],charge[i_pmt]); - } + for (int i_pmt=0;i_pmtFill(charge[detkey]); + if (time[detkey]!=0) time_PMTs->Fill(time[detkey]); + if (charge[detkey]!=0) charge_time_PMTs->Fill(time[detkey],charge[detkey]); + } //--------------------------------------------------------------- //------------- Determine max+min values ------------------------ //--------------------------------------------------------------- - maximum_pmts = 0; + maximum_pmts = 0; maximum_lappds = 0; maximum_time_pmts = 0; maximum_time_lappds = 0; min_time_pmts = 999.; min_time_lappds = 999.; - total_charge_pmts = 0; - for (int i_pmt=0;i_pmtmaximum_pmts) maximum_pmts = charge[i_pmt]; - total_charge_pmts+=charge[i_pmt]; - if (time[i_pmt]>maximum_time_pmts) maximum_time_pmts = time[i_pmt]; - if (time[i_pmt]maximum_pmts) maximum_pmts = charge[detkey]; + total_charge_pmts+=charge[detkey]; + if (time[detkey]>maximum_time_pmts) maximum_time_pmts = time[detkey]; + if (time[detkey]GetDetectorID(); - int LAPPDId = detectorkey_to_lappdid.at(detkey); std::vector& hits = apair.second; - if (verbose > 2) std::cout <<"detector key: "< 2) std::cout <<"detector key: "< temp_pos = ahit.GetPosition(); double x_lappd = temp_pos.at(0)-tank_center_x; @@ -506,16 +525,18 @@ bool EventDisplay::Execute(){ double z_lappd = temp_pos.at(2)-tank_center_z; double lappd_charge = 1.0; double t_lappd = ahit.GetTime(); - if ((lappds_selected && act_lappds[LAPPDId-1]==1) || !lappds_selected){ - if (verbose > 2) std::cout <<"Tube ID LAPPD: "< 2) std::cout <<"Detector Key LAPPD: "< 2) std::cout <<"LAPPD hit at x="<{}); } + time_lappd[detkey].push_back(t_lappd); if (t_lappd > maximum_time_lappds) maximum_time_lappds = t_lappd; if (t_lappd < min_time_lappds) min_time_lappds = t_lappd; } @@ -532,11 +553,19 @@ bool EventDisplay::Execute(){ //-------------------Fill LAPPD hists --------------------------- //--------------------------------------------------------------- - for (int i_lappd=0;i_lappdFill(charge_lappd[i_lappd]); - for (int i_time=0;i_timeFill(time_lappd[i_lappd].at(i_time)); + for (auto&& ahisto : charge_LAPPDs){ + unsigned long detkey = ahisto.first; + int histkey; + std::vector::iterator it = std::find(lappd_detkeys.begin(), lappd_detkeys.end(), detkey); + if (it != lappd_detkeys.end()) histkey = std::distance(lappd_detkeys.begin(), it); + else { + std::cout << "Detkey " << detkey <<" Not Found in Geometry class for LAPPDs!" << std::endl; + continue; + } + if ((lappds_selected && active_lappds_user.count(detkey)==1) || !lappds_selected){ + charge_LAPPDs[histkey]->Fill(charge_lappd[detkey]); + for (int i_time=0;i_timeFill(time_lappd[detkey].at(i_time)); } } } @@ -572,14 +601,14 @@ bool EventDisplay::Execute(){ leg_charge = new TLegend(0.6,0.6,0.9,0.9); leg_charge->SetLineColor(0); std::string lappd_str = "LAPPD "; - std::string lappd_nr = std::to_string(max_lappd_charge+1); + std::string lappd_nr = std::to_string(lappd_detkeys[max_lappd_charge]); std::string lappd_label = lappd_str+lappd_nr; leg_charge->AddEntry(charge_LAPPDs[max_lappd_charge],lappd_label.c_str()); for (int i_lappd=0;i_lappdDraw("same"); - std::string lappd_nr = std::to_string(i_lappd+1); + std::string lappd_nr = std::to_string(lappd_detkeys[i_lappd]); std::string lappd_label = lappd_str+lappd_nr; leg_charge->AddEntry(time_LAPPDs[i_lappd],lappd_label.c_str()); } @@ -591,14 +620,14 @@ bool EventDisplay::Execute(){ leg_time = new TLegend(0.6,0.6,0.9,0.9); //leg_time->SetNColumns(5); leg_time->SetLineColor(0); - lappd_nr = std::to_string(max_lappd_time+1); + lappd_nr = std::to_string(lappd_detkeys[max_lappd_time]); lappd_label = lappd_str+lappd_nr; leg_time->AddEntry(time_LAPPDs[max_lappd_time],lappd_label.c_str()); for (int i_lappd=1;i_lappdDraw("same"); - std::string lappd_nr = std::to_string(i_lappd+1); + std::string lappd_nr = std::to_string(lappd_detkeys[i_lappd]); std::string lappd_label = lappd_str+lappd_nr; leg_time->AddEntry(time_LAPPDs[i_lappd],lappd_label.c_str()); } @@ -1010,35 +1039,36 @@ void EventDisplay::draw_event_box(){ //calculate marker coordinates for (int i_pmt=0;i_pmtthreshold_time_high) || (threshold_time_low!=-999 && time[i_pmt]threshold_time_high) || (threshold_time_low!=-999 && time[detkey]SetMarkerColor(color_marker); marker_top->SetMarkerStyle(8); marker_top->SetMarkerSize(1); marker_pmts_top.push_back(marker_top); - } else if (fabs(y_pmt[i_pmt]-min_y)<0.01 || fabs(y_pmt[i_pmt]+1.30912)<0.01){ + } else if (fabs(y_pmt[detkey]-min_y)<0.01 || fabs(y_pmt[detkey]+1.30912)<0.01){ //draw PMTs on the bottom of tank - x[0]=0.5-size_top_drawing*x_pmt[i_pmt]/tank_radius; - y[0]=0.5-(tank_height/tank_radius+1)*size_top_drawing+size_top_drawing*z_pmt[i_pmt]/tank_radius; + x[0]=0.5-size_top_drawing*x_pmt[detkey]/tank_radius; + y[0]=0.5-(tank_height/tank_radius+1)*size_top_drawing+size_top_drawing*z_pmt[detkey]/tank_radius; TPolyMarker *marker_bottom = new TPolyMarker(1,x,y,""); - if (mode == "Charge") color_marker = Bird_Idx+int((charge[i_pmt]-threshold)/(maximum_pmts-threshold)*254); - else if (mode == "Time" && threshold_time_high == -999 && threshold_time_low == -999) color_marker = Bird_Idx+int((time[i_pmt]-min_time_overall)/(maximum_time_overall-min_time_overall)*254); - else if (mode == "Time" && threshold_time_low == -999 && threshold_time_high != -999) color_marker = Bird_Idx+int((time[i_pmt]-min_time_overall)/(threshold_time_high-min_time_overall)*254); - else if (mode == "Time" && threshold_time_low != -999 && threshold_time_high == -999) color_marker = Bird_Idx+int((time[i_pmt]-threshold_time_low)/(maximum_time_overall-threshold_time_low)*254); - else if (mode == "Time") color_marker = Bird_Idx+int((time[i_pmt]-threshold_time_low)/(threshold_time_high-threshold_time_low)*254); + if (mode == "Charge") color_marker = Bird_Idx+int((charge[detkey]-threshold)/(maximum_pmts-threshold)*254); + else if (mode == "Time" && threshold_time_high == -999 && threshold_time_low == -999) color_marker = Bird_Idx+int((time[detkey]-min_time_overall)/(maximum_time_overall-min_time_overall)*254); + else if (mode == "Time" && threshold_time_low == -999 && threshold_time_high != -999) color_marker = Bird_Idx+int((time[detkey]-min_time_overall)/(threshold_time_high-min_time_overall)*254); + else if (mode == "Time" && threshold_time_low != -999 && threshold_time_high == -999) color_marker = Bird_Idx+int((time[detkey]-threshold_time_low)/(maximum_time_overall-threshold_time_low)*254); + else if (mode == "Time") color_marker = Bird_Idx+int((time[detkey]-threshold_time_low)/(threshold_time_high-threshold_time_low)*254); //std::cout <<"time: "<SetMarkerColor(color_marker); marker_bottom->SetMarkerStyle(8); @@ -1048,23 +1078,23 @@ void EventDisplay::draw_event_box(){ //draw PMTs on the tank bulk double phi; - if (x_pmt[i_pmt]>0 && z_pmt[i_pmt]>0) phi = atan(z_pmt[i_pmt]/x_pmt[i_pmt])+TMath::Pi()/2; - else if (x_pmt[i_pmt]>0 && z_pmt[i_pmt]<0) phi = atan(x_pmt[i_pmt]/-z_pmt[i_pmt]); - else if (x_pmt[i_pmt]<0 && z_pmt[i_pmt]<0) phi = 3*TMath::Pi()/2+atan(z_pmt[i_pmt]/x_pmt[i_pmt]); - else if (x_pmt[i_pmt]<0 && z_pmt[i_pmt]>0) phi = TMath::Pi()+atan(-x_pmt[i_pmt]/z_pmt[i_pmt]); + if (x_pmt[detkey]>0 && z_pmt[detkey]>0) phi = atan(z_pmt[detkey]/x_pmt[detkey])+TMath::Pi()/2; + else if (x_pmt[detkey]>0 && z_pmt[detkey]<0) phi = atan(x_pmt[detkey]/-z_pmt[detkey]); + else if (x_pmt[detkey]<0 && z_pmt[detkey]<0) phi = 3*TMath::Pi()/2+atan(z_pmt[detkey]/x_pmt[detkey]); + else if (x_pmt[detkey]<0 && z_pmt[detkey]>0) phi = TMath::Pi()+atan(-x_pmt[detkey]/z_pmt[detkey]); else phi = 0.; if (phi>2*TMath::Pi()) phi-=(2*TMath::Pi()); phi-=TMath::Pi(); if (phi < - TMath::Pi()) phi = -TMath::Pi(); - if (phi<-TMath::Pi() || phi>TMath::Pi()) std::cout <<"Drawing Event: Phi out of bounds! "<<", x= "<SetMarkerColor(color_marker); marker_bulk->SetMarkerStyle(8); marker_bulk->SetMarkerSize(1); @@ -1095,14 +1125,17 @@ void EventDisplay::draw_event_box(){ double y[1]; double charge_single=1.; //FIXME when charge is implemented in LoadWCSimLAPPD bool draw_lappd_markers=false; - for (int i_lappd_hits=0;i_lappd_hits>::iterator lappd_hit_pos_it = hits_LAPPDs.begin(); + for(auto&& these_lappd_hit_positions : hits_LAPPDs){ + unsigned long detkey = these_lappd_hit_positions.first; + for (int i_single_lappd=0; i_single_lappdthreshold_time_high) || (threshold_time_low!=-999 && time_lappd_single0 && z_lappd>0) phi = atan(z_lappd/x_lappd)+TMath::Pi()/2; else if (x_lappd>0 && z_lappd<0) phi = atan(x_lappd/-z_lappd); else if (x_lappd<0 && z_lappd<0) phi = 3*TMath::Pi()/2+atan(z_lappd/x_lappd); diff --git a/UserTools/EventDisplay/EventDisplay.h b/UserTools/EventDisplay/EventDisplay.h index 97021f2f8..a193e9c79 100644 --- a/UserTools/EventDisplay/EventDisplay.h +++ b/UserTools/EventDisplay/EventDisplay.h @@ -104,7 +104,7 @@ class EventDisplay: public Tool { TBox *box_mrd_top = nullptr; double tank_height; double tank_radius; - static const int max_num_lappds = 200; //one is allowed to dream, right? + int max_num_lappds = 200; //one is allowed to dream, right? double detector_version; std::string detector_config; bool draw_ring_temp; @@ -116,8 +116,8 @@ class EventDisplay: public Tool { int n_mrd_pmts; int n_veto_pmts; int n_lappds; - std::vector active_lappds; - int act_lappds[max_num_lappds]={0}; + std::map active_lappds_user; // WCSim LAPPD IDs read from config file + std::map active_lappds; // converted to detectorkey double max_y; double min_y; @@ -143,18 +143,18 @@ class EventDisplay: public Tool { double tank_center_y; double tank_center_z; - std::map detectorkey_to_lappdid; + std::map lappdid_to_detectorkey; std::map channelkey_to_pmtid; std::map pmt_tubeid_to_channelkey; std::map channelkey_to_mrdpmtid; std::map channelkey_to_faccpmtid; static const unsigned long n_channels_per_lappd = 60; - double charge[200]; - double time[200]; - double charge_lappd[max_num_lappds]; - std::vector> hits_LAPPDs; - std::vector> time_lappd; + std::map charge; + std::map time; + std::map charge_lappd; + std::map> hits_LAPPDs; + std::map> time_lappd; std::vector mrddigittimesthisevent; std::vector mrddigitpmtsthisevent; std::vector mrddigitchargesthisevent; @@ -200,11 +200,14 @@ class EventDisplay: public Tool { int current_n_polylines = 0; //histograms + std::vector lappd_detkeys; + std::vector pmt_detkeys; + std::vector hitpmt_detkeys; TH1F *time_PMTs = nullptr; - TH1F *time_LAPPDs[max_num_lappds]; + std::map time_LAPPDs; //the keys for the LAPPD histogram maps are the indices of the lappd_detkeys vector TH1F *charge_PMTs = nullptr; TH2F *charge_time_PMTs = nullptr; - TH1F *charge_LAPPDs[max_num_lappds]; + std::map charge_LAPPDs; //legends TLegend *leg_charge = nullptr; diff --git a/UserTools/TankCalibrationDiffuser/.DS_Store b/UserTools/TankCalibrationDiffuser/.DS_Store deleted file mode 100644 index 5008ddfcf..000000000 Binary files a/UserTools/TankCalibrationDiffuser/.DS_Store and /dev/null differ diff --git a/UserTools/TankCalibrationDiffuser/README.md b/UserTools/TankCalibrationDiffuser/README.md index 4cf56f4a1..19b031fcc 100644 --- a/UserTools/TankCalibrationDiffuser/README.md +++ b/UserTools/TankCalibrationDiffuser/README.md @@ -14,16 +14,30 @@ The TankCalibrationDiffuser tool needs access to the hit information within ANNI The following variables can be configured for the TankCalibrationDiffuser tool: ``` -OutputFile /ANNIECode/ToolAnalysis/histogrammed_pmt_charge.root #Output root file for the current calibration run -StabilityFile /ANNIECode/ToolAnalysis/stability_pmtcharge.root #Output root file that displays the stability over multiple runs -GeometryFile configfiles/TankCalibrationDiffuser/geofile_128PMTs.txt #Geometry file specifying the different radii of the PMTs (taken from file geofile.txt in WCSim installation directory for the respective installation) -DiffuserX 0. #x-position of the diffuser ball -DiffuserY 0. #y-position of the diffuser ball -DiffuserZ 0. #z-position of the diffuser ball -ToleranceCharge 0.5 #tolerance of fit single p.e. value for being classified as a bad PMT -ToleranceTime 0.5 #tolerance of mean time value [ns] for being classified as a bad PMT -TApplication 0 #0/1, depending on whether plots should be shown interactively or not +# TankCalibrationDiffuser Config File + +OutputFile SimulationCalibrationTest #Output root prefix name for the current run +DiffuserX 0. #x-position of the diffuser ball +DiffuserY 0. #y-position of the diffuser ball +DiffuserZ 0. #z-position of the diffuser ball +ToleranceCharge 0.5 #tolerance of fit single p.e. value for being classified as a bad PMT +ToleranceTime 0.5 #tolerance of mean time value [ns] for being classified as a bad PMT +FitMethod Gaus2Exp #fit function for charge, options: Gaus2Exp (2 times gaus + exp), Gaus2 (2 times gaus), Gaus (single gaus) +TApplication 0 #0/1, depending on whether plots should be shown interactively or not verbose 1 #verbosity of the application ``` + +## OutputFiles + +The tool produces two output files: +* a root file `_PMTStability_Run.root` which contains the charge & time histograms for all the PMTs and overview plots of the fitted charge and time distributions +* a txt file `_Run_pmts_laser_calibration.txt`, which contains a PMT-by-PMT summary of the calibration run information. The columns are (from left to right) + * PMT Detkey + * Fitted charge - mean + * Fitted charge - sigma + * Fitted time - mean + * Fitted time - sigma + * Time deviation (observed / expected hit time) +* The last line of the .txt-file contains the average fit information averaged over all PMTs (assigned to detectorkey 10000) diff --git a/UserTools/TankCalibrationDiffuser/TankCalibrationDiffuser.cpp b/UserTools/TankCalibrationDiffuser/TankCalibrationDiffuser.cpp index 5a242f579..42bf8ba24 100644 --- a/UserTools/TankCalibrationDiffuser/TankCalibrationDiffuser.cpp +++ b/UserTools/TankCalibrationDiffuser/TankCalibrationDiffuser.cpp @@ -17,23 +17,22 @@ bool TankCalibrationDiffuser::Initialise(std::string configfile, DataModel &data if(configfile!="") m_variables.Initialise(configfile); //loading config file m_data= &data; //assigning transient data pointer - + //---------------------------------------------------------------------------- //---------------Get configuration variables for this tool-------------------- //---------------------------------------------------------------------------- - + m_variables.Get("OutputFile",outputfile); - m_variables.Get("StabilityFile",stabilityfile); - m_variables.Get("GeometryFile",geometryfile); m_variables.Get("DiffuserX",diffuser_x); m_variables.Get("DiffuserY",diffuser_y); m_variables.Get("DiffuserZ",diffuser_z); m_variables.Get("ToleranceCharge",tolerance_charge); m_variables.Get("ToleranceTime",tolerance_time); + m_variables.Get("FitMethod",FitMethod); m_variables.Get("TApplication",use_tapplication); m_variables.Get("verbose",verbose); - help_file = new TFile("help_file.root","RECREATE"); //let the histos be associated with this file + help_file = new TFile("configfiles/TankCalibrationDiffuser/help_file.root","RECREATE"); //let the histograms be associated with this file //---------------------------------------------------------------------------- //---------------Get basic geometry properties ------------------------------- @@ -57,45 +56,28 @@ bool TankCalibrationDiffuser::Initialise(std::string configfile, DataModel &data //---------------Read in geometry properties of PMTs-------------------------- //---------------------------------------------------------------------------- - if (verbose > 0) std::cout <<"Geofile name: "<::iterator it = Detectors->at("Tank").begin(); - it != Detectors->at("Tank").end(); - ++it){ + for (std::map::iterator it = Detectors->at("Tank").begin(); it != Detectors->at("Tank").end(); ++it){ Detector* apmt = it->second; unsigned long detkey = it->first; unsigned long chankey = apmt->GetChannels()->begin()->first; + pmt_detkeys.push_back(detkey); + std::string dettype = apmt->GetDetectorType(); + + if (verbose > 1) std::cout <<"PMT with detkey: "< 1) std::cout <<"PMT Print:" <Print()<GetDetectorPosition(); if (verbose > 2) std::cout <<"detkey: "<(detkey,position_PMT.Y()-tank_center_y)); z_PMT.insert(std::pair(detkey,position_PMT.Z()-tank_center_z)); - double rho = sqrt(pow(x_PMT.at(detkey),2)+pow(y_PMT.at(detkey),2)); if (x_PMT.at(detkey) < 0) rho*=-1; rho_PMT.insert(std::pair(detkey,rho)); - - double phi; if (x_PMT.at(detkey)>0 && z_PMT.at(detkey)>0) phi = atan(x_PMT.at(detkey)/z_PMT.at(detkey)); if (x_PMT.at(detkey)>0 && z_PMT.at(detkey)<0) phi = TMath::Pi()/2+atan(x_PMT.at(detkey)/-z_PMT.at(detkey)); @@ -118,51 +97,63 @@ bool TankCalibrationDiffuser::Initialise(std::string configfile, DataModel &data if (x_PMT.at(detkey)<0 && z_PMT.at(detkey)>0) phi = 3*TMath::Pi()/2+atan(-x_PMT.at(detkey)/z_PMT.at(detkey)); phi_PMT.insert(std::pair(detkey,phi)); - - double expectedT = (sqrt(pow(x_PMT.at(detkey)-diffuser_x,2)+pow(y_PMT.at(detkey)-diffuser_y,2)+pow(z_PMT.at(detkey)-diffuser_z,2))-radius_PMT[detkey])/c_vacuum*n_water*1E9; - expected_time.insert(std::pair(detkey,expectedT)); if (verbose > 2) std::cout <<"detectorkey: "< 2) std::cout <<"rho PMT "< 2) std::cout <<"y PMT: "<(detkey,0)); + + if (dettype.find("R5912") != std::string::npos) radius_PMT[detkey] = map_type_radius["R5912HQE"]; + else if (dettype.find("R7081") != std::string::npos) radius_PMT[detkey] = map_type_radius["R7081"]; + else if (dettype.find("D784KFLB") != std::string::npos) radius_PMT[detkey] = map_type_radius["D784KFLB"]; + else if (dettype.find("EMI9954KB") != std::string::npos) radius_PMT[detkey] = map_type_radius["EMI9954KB"]; + + double expectedT = (sqrt(pow(x_PMT.at(detkey)-diffuser_x,2)+pow(y_PMT.at(detkey)-diffuser_y,2)+pow(z_PMT.at(detkey)-diffuser_z,2))-radius_PMT[detkey])/c_vacuum*n_water*1E9; + expected_time.insert(std::pair(detkey,expectedT)); } - std::cout <<"Number of tank PMTs: "< 1) std::cout <<"Number of tank PMTs: "< 1) std::cout <<"Detkey: "< charge_hist; hist_charge = new TH1F("hist_charge","Overall charge distribution (all PMTs)",500,1,0); hist_time = new TH1F("hist_time","Overall time distribution (all PMTs)",500,1,0); hist_tubeid = new TH1F("hist_tubeid","Overall Tube ID distribution",500,1,0); hist_charge_2D_y_phi = new TH2F("hist_charge_2D_y_phi","Spatial distribution of charge (all PMTs)",100,0,360,25,-2.5,2.5); - hist_time_2D_y_phi = new TH2F("hist_time_2D_y_phi","Spatial distribution of time (all PMTs)",100,0,360,25,-2.5,2.5); + hist_time_2D_y_phi = new TH2F("hist_time_2D_y_phi","Spatial distribution of time deviations (all PMTs)",100,0,360,25,-2.5,2.5); + hist_time_2D_y_phi_mean = new TH2F("hist_time_2D_y_phi_mean","Spatial distribution of time (all PMTs)",100,0,360,25,-2.5,2.5); for (int i_tube=0;i_tubecd(); //---------------------------------------------------------------------------- //---Make histograms appear during execution by declaring TApplication-------- //---------------------------------------------------------------------------- - + if (use_tapplication){ - if (verbose > 0) std::cout <<"open TApplication..."< 0) std::cout <<"TApplication initialized..."< 0) std::cout <<"Opening TApplication..."< 0) std::cout <<"TApplication initialized..."< 0) std::cout <<"Executing Tool TankCalibrationDiffuser ..."<Stores.count("ANNIEEvent"); if(!annieeventexists){ std::cerr<<"Error: No ANNIEEvent store!"<Stores["ANNIEEvent"]->Get("EventNumber", evnum); m_data->Stores["ANNIEEvent"]->Get("BeamStatus",BeamStatus); - //---------------------------------------------------------------------------- //---------------Read out charge and hit values of PMTs----------------------- //---------------------------------------------------------------------------- - + for (int i_pmt = 0; i_pmt < n_tank_pmts; i_pmt++){ - PMT_ishit.at(i_pmt) = 0; + unsigned long detkey = pmt_detkeys[i_pmt]; + PMT_ishit[detkey] = 0; } int vectsize = MCHits->size(); @@ -230,22 +221,22 @@ bool TankCalibrationDiffuser::Execute(){ int detectorkey = thistube->GetDetectorID(); if (thistube->GetDetectorElement()=="Tank"){ std::vector& Hits = apair.second; - int wcsim_id; PMT_ishit[detectorkey] = 1; for (MCHit &ahit : Hits){ - if (verbose > 2) std::cout <<"charge "< 2) std::cout <<"Charge "<Fill(ahit.GetCharge()); hist_time->Fill(ahit.GetTime()); - hist_tubeid->Fill(ahit.GetTubeId()); - hist_charge_singletube[ahit.GetTubeId()]->Fill(ahit.GetCharge()); - hist_time_singletube[ahit.GetTubeId()]->Fill(ahit.GetTime()); - } + hist_tubeid->Fill(detectorkey); + hist_charge_singletube[detectorkey]->Fill(ahit.GetCharge()); + hist_time_singletube[detectorkey]->Fill(ahit.GetTime()); + } } } //check whether PMT_ishit is filled correctly for (int i_pmt = 0; i_pmt < n_tank_pmts ; i_pmt++){ - if (verbose > 2) std::cout <<"PMT "< 2) std::cout <<"PMT "< 0) std::cout <<"TankCalibrationDiffuser: Finalise"<Print(); + //hist_charge->Print(); hist_charge->Write(); hist_time->Write(); hist_tubeid->Write(); // create txt file to store calibration data for the specific run - std::string filename_pre = "Run"; + std::string filename_pre = "_Run"; std::string filename_post = "_pmts_laser_calibration.txt"; - std::string filename_complete = filename_pre+newRunNumber+filename_post; + std::string filename_complete = outputfile+filename_pre+newRunNumber+filename_post; std::cout <<"writing to txt file: "<GetMean(); - Double_t par[8]={10,0.3,0.1,10,1.0,0.5,1,-1}; - TF1 *total = new TF1("total","gaus(0)+gaus(3)+expo(6)",0,10); + unsigned long detkey = pmt_detkeys[i_tube]; + double mean_charge=hist_charge_singletube[detkey]->GetMean(); + + Double_t par_gaus2exp[8] = {10,0.3,0.1,10,1.0,0.5,1,-1}; + Double_t par_gaus2[6] = {10,0.3,0.1,10,1.0,0.5}; + Double_t par_gaus[3] = {10,1.0,0.5}; + + TF1 *total; + if (FitMethod == "Gaus2Exp") total = new TF1("total","gaus(0)+gaus(3)+expo(6)",0,10); + else if (FitMethod == "Gaus2") total = new TF1("total","gaus(0)+gaus(3)",0,10); + else if (FitMethod == "Gaus") total = new TF1("total","gaus",0,10); + else { + std::cout <<"ERROR (TankCalibrationDiffuser): FitFunction is not part of the options, please extend the options Using standard Gaus."<SetLineColor(2); - - total->SetParameters(par); - + if (FitMethod == "Gaus2Exp") total->SetParameters(par_gaus2exp); + else if (FitMethod == "Gaus2") total->SetParameters(par_gaus2); + else total->SetParameters(par_gaus2); + //total = new TF1("total",total_value,0,10,8); //total->SetParameters(10,0.3,0.1,10,1.,0.3,1,1); - hist_charge_singletube[i_tube]->Fit(total,"QR+"); + TFitResultPtr FitResult = hist_charge_singletube[detkey]->Fit(total,"QR+"); + Int_t FitResult_int = FitResult; //hist_charge_singletube[i_tube]->Fit("gaus[0]+gaus[3]+exp[6]"); - hist_charge_singletube[i_tube]->Write(); - TF1 *fit_result_charge=hist_charge_singletube[i_tube]->GetFunction("total"); - mean_charge_fit[i_tube]=fit_result_charge->GetParameter(4); - hist_charge_mean->Fill(mean_charge_fit[i_tube]); - - - double mean_time=hist_time_singletube[i_tube]->GetMean(); - hist_time_singletube[i_tube]->Fit("gaus","Q"); - hist_time_singletube[i_tube]->Write(); - TF1 *fit_result_time=hist_time_singletube[i_tube]->GetFunction("gaus"); - mean_time_fit[i_tube]=fit_result_time->GetParameter(1); - hist_time_mean->Fill(mean_time_fit[i_tube]-expected_time[i_tube]); + hist_charge_singletube[detkey]->Write(); + if (FitResult_int ==0){ + TF1 *fit_result_charge=hist_charge_singletube[detkey]->GetFunction("total"); + if (FitMethod == "Gaus2Exp") { + mean_charge_fit[detkey] = fit_result_charge->GetParameter(4); + rms_charge_fit[detkey] = fit_result_charge->GetParameter(5); + } + else if (FitMethod == "Gaus2") { + mean_charge_fit[detkey]=fit_result_charge->GetParameter(4); + mean_charge_fit[detkey]=fit_result_charge->GetParameter(4); + } + else { + mean_charge_fit[detkey]=fit_result_charge->GetParameter(1); + rms_charge_fit[detkey]=fit_result_charge->GetParameter(2); + } + } else { + mean_charge_fit[detkey] = 0.; + rms_charge_fit[detkey] = 0.; + } + hist_charge_mean->Fill(mean_charge_fit[detkey]); + + double mean_time=hist_time_singletube[detkey]->GetMean(); + TFitResultPtr FitResultTime = hist_time_singletube[detkey]->Fit("gaus","Q"); //time fit is always the same for now, assume simple Gaussian + Int_t FitResultTime_int = FitResultTime; + hist_time_singletube[detkey]->Write(); + if (FitResultTime_int == 0) { //fit was okay and has a result + TF1 *fit_result_time=hist_time_singletube[detkey]->GetFunction("gaus"); + mean_time_fit[detkey]=fit_result_time->GetParameter(1); + rms_time_fit[detkey]=fit_result_time->GetParameter(2); + } + else { + mean_time_fit[detkey] = 0.; + rms_time_fit[detkey] = 0.; + } + hist_time_dev->Fill(mean_time_fit[detkey]-expected_time[detkey]); + hist_time_mean->Fill(mean_time_fit[detkey]); if (verbose > 2){ - std::cout <<"expected hit time: "< set to 0 + mean_charge_fit[detkey] = 0.; + rms_charge_fit[detkey] = 0.; + } - if (hist_time_2D_y_phi->GetBinContent(int(phi_PMT[i_tube]/2/TMath::Pi()*100)+1,int((y_PMT[i_tube]+2.5)/5.*25)+1)==0){ - hist_time_2D_y_phi->SetBinContent(int(phi_PMT[i_tube]/2/TMath::Pi()*100)+1,int((y_PMT[i_tube]+2.5)/5.*25)+1,fabs(mean_time_fit[i_tube]-expected_time[i_tube])); - hist_charge_2D_y_phi->SetBinContent(int(phi_PMT[i_tube]/2/TMath::Pi()*100)+1,int((y_PMT[i_tube]+2.5)/5.*25)+1,mean_charge_fit[i_tube]); - } - else { - hist_time_2D_y_phi->SetBinContent(int(phi_PMT[i_tube]/2/TMath::Pi()*100)+1,int((y_PMT[i_tube]+2.5)/5.*25)+2,fabs(mean_time_fit[i_tube]-expected_time[i_tube])); - hist_charge_2D_y_phi->SetBinContent(int(phi_PMT[i_tube]/2/TMath::Pi()*100)+1,int((y_PMT[i_tube]+2.5)/5.*25)+2,mean_charge_fit[i_tube]); - } - result_file<GetMean()<<" "<GetRMS()<<" "<GetMean()<<" "<GetRMS()<GetBinContent(int(phi_PMT[detkey]/2/TMath::Pi()*100)+1,int((y_PMT[detkey]+2.5)/5.*25)+1)==0){ + hist_time_2D_y_phi->SetBinContent(int(phi_PMT[detkey]/2/TMath::Pi()*100)+1,int((y_PMT[detkey]+2.5)/5.*25)+1,fabs(mean_time_fit[detkey]-expected_time[detkey])); + hist_time_2D_y_phi_mean->SetBinContent(int(phi_PMT[detkey]/2/TMath::Pi()*100)+1,int((y_PMT[detkey]+2.5)/5.*25)+1,mean_time_fit[detkey]); + hist_charge_2D_y_phi->SetBinContent(int(phi_PMT[detkey]/2/TMath::Pi()*100)+1,int((y_PMT[detkey]+2.5)/5.*25)+1,mean_charge_fit[detkey]); + } + else { + hist_time_2D_y_phi->SetBinContent(int(phi_PMT[detkey]/2/TMath::Pi()*100)+1,int((y_PMT[detkey]+2.5)/5.*25)+2,fabs(mean_time_fit[detkey]-expected_time[detkey])); + hist_time_2D_y_phi_mean->SetBinContent(int(phi_PMT[detkey]/2/TMath::Pi()*100)+1,int((y_PMT[detkey]+2.5)/5.*25)+2,mean_time_fit[detkey]); + hist_charge_2D_y_phi->SetBinContent(int(phi_PMT[detkey]/2/TMath::Pi()*100)+1,int((y_PMT[detkey]+2.5)/5.*25)+2,mean_charge_fit[detkey]); + } + result_file<GetMean()<<" "<GetRMS()<<" "<GetMean()<<" "<GetRMS()<<" "<GetMean()<<" "<GetRMS()<GetXaxis()->SetTitle("t_{arrival} - t_{expected} [ns]"); + hist_time_mean->GetXaxis()->SetTitle("t_{arrival} [ns]"); + hist_time_dev->GetXaxis()->SetTitle("t_{arrival} - t_{expected} [ns]"); hist_charge_mean->GetXaxis()->SetTitle("charge [p.e.]"); hist_charge_2D_y_phi->GetXaxis()->SetTitle("#phi [deg]"); hist_charge_2D_y_phi->GetYaxis()->SetTitle("y [m]"); + hist_charge_2D_y_phi->SetStats(0); hist_time_2D_y_phi->GetXaxis()->SetTitle("#phi [deg]"); hist_time_2D_y_phi->GetYaxis()->SetTitle("y [m]"); hist_time_2D_y_phi->SetStats(0); - hist_charge_2D_y_phi->SetStats(0); + hist_time_2D_y_phi_mean->GetXaxis()->SetTitle("#phi [deg]"); + hist_time_2D_y_phi_mean->GetYaxis()->SetTitle("y [m]"); + hist_time_2D_y_phi_mean->SetStats(0); hist_time_mean->Write(); + hist_time_dev->Write(); hist_charge_mean->Write(); hist_time_2D_y_phi->Write(); + hist_time_2D_y_phi_mean->Write(); hist_charge_2D_y_phi->Write(); - + //---------------------------------------------------------------------------- //---------------Create and write stability plots------------------------ //---------------------------------------------------------------------------- @@ -359,56 +402,69 @@ bool TankCalibrationDiffuser::Finalise(){ gr_stability_time = new TGraphErrors(); gr_stability_time->SetName("gr_stability_time"); gr_stability_time->SetTitle("Stability PMT time calibration"); + gr_stability_time_mean = new TGraphErrors(); + gr_stability_time_mean->SetName("gr_stability_time_mean"); + gr_stability_time_mean->SetTitle("Stability PMT time (mean)"); + + // + //read in information from 100 last runs to produce stability/time evolution plots (if they do not exist, just skip them...) + // - //read in information from 100 last runs (if they do not exist, just skip them...) const int n_entries=100; double run_numbers[n_entries]; double entries_charge[n_entries]; double rms_charge[n_entries]; double entries_time[n_entries]; double rms_time[n_entries]; + double entries_time_mean[n_entries]; + double rms_time_mean[n_entries]; bool file_exists[n_entries]; int runnumber_temp; - + // //get stability values for time and charge from txt-files + // + for (int i_run=n_entries-1;i_run>=0;i_run--){ - + std::stringstream ss_runnumber; runnumber_temp = runnumber-i_run; ss_runnumber< 1) std::cout <<"Stability: Filename: "<>dummy_temp>>dummy_temp>>dummy_temp; - else file_temp>>dummy_temp>>entries_charge[n_entries-1-i_run]>>rms_charge[n_entries-1-i_run]>>entries_time[n_entries-1-i_run]>>rms_time[n_entries-1-i_run]; + if (i_pmt>dummy_temp>>dummy_temp>>dummy_temp>>dummy_temp>>dummy_temp>>dummy_temp; + else file_temp>>dummy_temp>>entries_charge[n_entries-1-i_run]>>rms_charge[n_entries-1-i_run]>>entries_time_mean[n_entries-1-i_run]>>rms_time_mean[n_entries-1-i_run]>>entries_time[n_entries-1-i_run]>>rms_time[n_entries-1-i_run]; } file_temp.close(); } else { file_exists[n_entries-1-i_run]=false; } - } - int i_point_graph=0; + // //fill TGraphErrors with stability values for time and charge + // + + int i_point_graph=0; for (int i_run=n_entries-1;i_run>=0;i_run--){ if (file_exists[n_entries-1-i_run]){ - if (verbose > 1) std::cout <<"Stability, run "< 1) std::cout <<"Stability, Run "<SetPoint(i_point_graph,runnumber-i_run,entries_charge[n_entries-1-i_run]); gr_stability->SetPointError(i_point_graph,0.,rms_charge[n_entries-1-i_run]); gr_stability_time->SetPoint(i_point_graph,runnumber-i_run,entries_time[n_entries-1-i_run]); gr_stability_time->SetPointError(i_point_graph,0.,rms_time[n_entries-1-i_run]); + gr_stability_time_mean->SetPoint(i_point_graph,runnumber-i_run,entries_time_mean[n_entries-1-i_run]); + gr_stability_time_mean->SetPointError(i_point_graph,0.,rms_time_mean[n_entries-1-i_run]); i_point_graph++; } - } gr_stability->SetMarkerStyle(22); @@ -423,10 +479,20 @@ bool TankCalibrationDiffuser::Finalise(){ gr_stability_time->SetLineWidth(2); gr_stability_time->GetXaxis()->SetTitle("Run Number"); gr_stability_time->GetYaxis()->SetTitle("hit time deviaton [ns]"); + gr_stability_time_mean->SetMarkerStyle(22); + gr_stability_time_mean->SetMarkerSize(1.0); + gr_stability_time_mean->SetLineColor(1); + gr_stability_time_mean->SetLineWidth(2); + gr_stability_time_mean->GetXaxis()->SetTitle("Run Number"); + gr_stability_time_mean->GetYaxis()->SetTitle("hit time [ns]"); gr_stability->Write("gr_stability"); gr_stability_time->Write("gr_stability_time"); + gr_stability_time_mean->Write("gr_stability_time_mean"); + + // + //Fill new entries into histogram + // - //fill new entries into histogram std::map::iterator it_time = bad_time.begin(); std::map::iterator it_charge = bad_charge.begin(); @@ -437,146 +503,167 @@ bool TankCalibrationDiffuser::Finalise(){ std::cout <<"///////////////////////////////////////////////////////////////////"<tolerance_time) { + std::cout <<"Abnormally high time deviation for detkey "<(detkey,mean_time_fit[detkey]-expected_time[detkey])); + } + if (mean_charge_fit[detkey]<(1.0-tolerance_charge) || mean_charge_fit[detkey]>(1.0+tolerance_charge)){ + std::cout <<"Abnormal charge for detkey "<(detkey,mean_charge_fit[detkey])); + } + } + + //---------------------------------------------------------------------------- + //---------------Save stability root plots------------------------------------ + //---------------------------------------------------------------------------- - if (fabs(mean_time_fit[i_tube]-expected_time[i_tube])>tolerance_time) { - std::cout <<"Abnormally high time deviation for tube "<(i_tube,mean_time_fit[i_tube]-expected_time[i_tube])); + canvas_overview = new TCanvas("canvas_overview","Stability PMTs (Deviations)",900,600); + canvas_overview->Divide(3,2); + canvas_overview->cd(1); + hist_charge_2D_y_phi->Draw("colz"); + std::map::iterator it_charge_read; //draw red boxes around the misaligned PMTs (in charge as well as in time) + if (verbose > 0) std::cout << "PMTs with bad charge: "<first; + int bin_charge_x = hist_charge_2D_y_phi->GetXaxis()->FindBin(phi_PMT[pmt]/TMath::Pi()*180.); + int bin_charge_y = hist_charge_2D_y_phi->GetYaxis()->FindBin(y_PMT[pmt]); + if (verbose > 0) std::cout << "PMT: "<(1.0+tolerance_charge)){ - std::cout <<"Abnormal charge for tube "<(i_tube,mean_charge_fit[i_tube])); + canvas_overview->cd(2); + hist_charge_mean->Draw(); + canvas_overview->cd(3); + gr_stability->Draw("AEL"); + canvas_overview->cd(4); + hist_time_2D_y_phi->Draw("colz"); + std::map::iterator it_time_read; + if (verbose > 0) std::cout << "PMTs with bad timing: "<first; + //if (verbose > 0) std::cout <<"PMT: "<GetXaxis()->FindBin(phi_PMT[pmt]/TMath::Pi()*180.); + int bin_time_y = hist_time_2D_y_phi->GetYaxis()->FindBin(y_PMT[pmt]); + if (verbose > 0) std::cout << "PMT: "<cd(5); + hist_time_dev->Draw(); + canvas_overview->cd(6); + gr_stability_time->Draw("ALE"); + canvas_overview->Modified(); + canvas_overview->Update(); + canvas_overview->Write(); + //---------------------------------------------------------------------------- - //---------------Save stability root plots------------------------------------ + //---------------Save stability root plots v2 (no deviations)----------------- //---------------------------------------------------------------------------- - - canvas_overview = new TCanvas("canvas_overview","Stability PMTs",900,600); - canvas_overview->Divide(3,2); - canvas_overview->cd(1); - hist_charge_2D_y_phi->Draw("colz"); - std::map::iterator it_charge_read; //draw red boxes around the misaligned PMTs (in charge as well as in time) - if (verbose > 0) std::cout << "PMTs with bad charge: "<first; - int bin_charge_x = hist_charge_2D_y_phi->GetXaxis()->FindBin(phi_PMT[pmt]/TMath::Pi()*180.); - int bin_charge_y = hist_charge_2D_y_phi->GetYaxis()->FindBin(y_PMT[pmt]); - if (verbose > 0) std::cout << "PMT: "<cd(2); - hist_charge_mean->Draw(); - canvas_overview->cd(3); - gr_stability->Draw("AEL"); - canvas_overview->cd(4); - hist_time_2D_y_phi->Draw("colz"); - std::map::iterator it_time_read; - if (verbose > 0) std::cout << "PMTs with bad timing: "<first; - //if (verbose > 0) std::cout <<"PMT: "<GetXaxis()->FindBin(phi_PMT[pmt]/TMath::Pi()*180.); - int bin_time_y = hist_time_2D_y_phi->GetYaxis()->FindBin(y_PMT[pmt]); - if (verbose > 0) std::cout << "PMT: "<cd(5); - hist_time_mean->Draw(); - canvas_overview->cd(6); - gr_stability_time->Draw("ALE"); - canvas_overview->Modified(); - canvas_overview->Update(); - canvas_overview->Write(); - - file_out->Close(); + + canvas_overview2 = new TCanvas("canvas_overview2","Stability PMTs (Mean)",900,600); + canvas_overview2->Divide(3,2); + canvas_overview2->cd(1); + hist_charge_2D_y_phi->Draw("colz"); + canvas_overview2->cd(2); + hist_charge_mean->Draw(); + canvas_overview2->cd(3); + gr_stability->Draw("AEL"); + canvas_overview2->cd(4); + hist_time_2D_y_phi_mean->Draw("colz"); + canvas_overview2->cd(5); + hist_time_mean->Draw(); + canvas_overview2->cd(6); + gr_stability_time_mean->Draw("ALE"); + canvas_overview2->Modified(); + canvas_overview2->Update(); + canvas_overview2->Write(); + + file_out->Close(); //---------------------------------------------------------------------------- //---------------Prompt user inputs (to avoid abortion of canvas)------------- //---------------------------------------------------------------------------- - int valid_calibration; - if (verbose > 0) { - std::cout << "//////////////////////////////////////////////////////////////"<>valid_calibration; - app_stability->Terminate(); - } - } - - help_file->Close(); + int valid_calibration; + if (verbose > 0) { + std::cout << "//////////////////////////////////////////////////////////////"<>valid_calibration; + app_stability->Terminate(); + } + } + + help_file->Close(); //---------------------------------------------------------------------------- //------------------------Deleting all objects ------------------------------- //---------------------------------------------------------------------------- - //histograms deleted by closing the files ealier - - std::cout <<"delete TBoxes..."< 0) std::cout <<"Deleting TBoxes..."< 0) std::cout <<"Deleting TF1s..."< 0) std::cout <<"Deleting canvas, application, files, geom, mchits"<Print();*/ - return true; + return true; } - //---------------------------------------------------------------------------- - //---------------------Helper functions -------------------------------------- - //---------------------------------------------------------------------------- - - bool TankCalibrationDiffuser::does_file_exist(const char *fileName){ - std::ifstream infile(fileName); - return infile.good(); - } - - void TankCalibrationDiffuser::draw_red_box(TH2 *hist, int bin_x, int bin_y, int verbosity){ +//---------------------------------------------------------------------------- +//---------------------Helper functions -------------------------------------- +//---------------------------------------------------------------------------- - Int_t bin_global = hist->GetBin(bin_x, bin_y,0); +bool TankCalibrationDiffuser::does_file_exist(const char *fileName){ + std::ifstream infile(fileName); + return infile.good(); +} - double size_x=2.0; - double size_y=0.05; - TBox *b = new TBox(hist->GetXaxis()->GetBinLowEdge(bin_x)-size_x,hist->GetYaxis()->GetBinLowEdge(bin_y)-size_y,hist->GetXaxis()->GetBinWidth(bin_x)+hist->GetXaxis()->GetBinLowEdge(bin_x)+size_x, hist->GetYaxis()->GetBinWidth(bin_y)+hist->GetYaxis()->GetBinLowEdge(bin_y)+size_y); - - if (verbosity > 2) { - std::cout <<"Draw Red Box: lower x: "<GetXaxis()->GetBinLowEdge(bin_x)-0.1<GetYaxis()->GetBinLowEdge(bin_y)-0.1<GetXaxis()->GetBinWidth(bin_x)+hist->GetXaxis()->GetBinLowEdge(bin_x)+0.1<GetXaxis()->GetBinWidth(bin_y)+hist->GetXaxis()->GetBinLowEdge(bin_y)+0.1<SetFillStyle(0); - b->SetLineColor(2); - b->SetLineWidth(2); - b->Draw(); + Int_t bin_global = hist->GetBin(bin_x, bin_y,0); - vector_tbox.push_back(b); + double size_x=2.0; + double size_y=0.05; + TBox *b = new TBox(hist->GetXaxis()->GetBinLowEdge(bin_x)-size_x,hist->GetYaxis()->GetBinLowEdge(bin_y)-size_y,hist->GetXaxis()->GetBinWidth(bin_x)+hist->GetXaxis()->GetBinLowEdge(bin_x)+size_x, hist->GetYaxis()->GetBinWidth(bin_y)+hist->GetYaxis()->GetBinLowEdge(bin_y)+size_y); + + if (verbosity > 2) { + std::cout <<"Draw Red Box: lower x: "<GetXaxis()->GetBinLowEdge(bin_x)-0.1<GetYaxis()->GetBinLowEdge(bin_y)-0.1<GetXaxis()->GetBinWidth(bin_x)+hist->GetXaxis()->GetBinLowEdge(bin_x)+0.1<GetXaxis()->GetBinWidth(bin_y)+hist->GetXaxis()->GetBinLowEdge(bin_y)+0.1<SetFillStyle(0); + b->SetLineColor(2); + b->SetLineWidth(2); + b->Draw(); - + vector_tbox.push_back(b); +} diff --git a/UserTools/TankCalibrationDiffuser/TankCalibrationDiffuser.h b/UserTools/TankCalibrationDiffuser/TankCalibrationDiffuser.h index 8daf1b1f3..1cf0c7679 100644 --- a/UserTools/TankCalibrationDiffuser/TankCalibrationDiffuser.h +++ b/UserTools/TankCalibrationDiffuser/TankCalibrationDiffuser.h @@ -27,97 +27,110 @@ #include "TH2F.h" #include "TFile.h" +/** + * \class MonitorTankLive +* +* $Author: M. Nieslony $ +* $Date: 2019/08/09 12:55:00 $ +* Contact: mnieslon@uni-mainz.de +*/ + class TankCalibrationDiffuser: public Tool { - public: - - TankCalibrationDiffuser(); - bool Initialise(std::string configfile,DataModel &data); - bool Execute(); - bool Finalise(); - - bool does_file_exist(const char *fileName); - void draw_red_box(TH2 *hist, int bin_x, int bin_y, int verbosity); - - private: - - //define input variables - std::string outputfile; - std::string stabilityfile; - std::string geometryfile; - double diffuser_x; - double diffuser_y; - double diffuser_z; - double tolerance_charge; - double tolerance_time; - bool use_tapplication; - int verbose; - - //define ANNIEEvent variables - int evnum; - int runnumber; - TimeClass* EventTime=nullptr; - std::vector* TriggerData; - BeamStatusClass* BeamStatus=nullptr; - std::map>* MCHits = nullptr; - Geometry *geom = nullptr; - Detector det; - std::map channelkey_to_pmtid; - std::map pmt_tubeid_to_channelkey; - - - //define useful variables - const double n_water = 1.33; - const double c_vacuum= 2.99792E8; //in m/s - int n_tank_pmts, n_mrd_pmts, n_veto_pmts, n_lappds; - double tank_center_x, tank_center_y, tank_center_z; - static const unsigned long n_channels_per_lappd = 60; - - //define monitoring variables - std::map bad_time, bad_charge; - - //define arrays for storing PMT calibration values - std::map x_PMT; - std::map y_PMT; - std::map z_PMT; - std::map rho_PMT; - std::map phi_PMT; - std::map radius_PMT; - std::map PMT_ishit; - std::map mean_charge_fit; - std::map mean_time_fit; - std::map expected_time; - - //define histograms for stability plots - TH1F *hist_tubeid = nullptr; - TH1F *hist_charge = nullptr; - TH1F *hist_time = nullptr; - std::map hist_charge_singletube; - std::map hist_time_singletube; - TH1F *hist_charge_mean = nullptr; - TH1F *hist_time_mean = nullptr; - TH2F *hist_charge_2D_y_phi = nullptr; - TH2F *hist_time_2D_y_phi = nullptr; - - //container for red boxes surrounding badly calibrated PMTs - std::vector vector_tbox; - - //container for TF1s used to fit the distributions - std::vector vector_tf1; - - //define graphs for stability plots - TGraphErrors *gr_stability = nullptr; - TGraphErrors *gr_stability_time = nullptr; - - //define file to save data - TFile *file_out = nullptr; - TFile *help_file = nullptr; - - //define overview canvas - TCanvas *canvas_overview = nullptr; - - //define TApplication for possibility of interactive plotting - TApplication *app_stability = nullptr; + public: + + TankCalibrationDiffuser(); ///< Simple constructor + bool Initialise(std::string configfile,DataModel &data); ///< Initialise Function for setting up Tool resources. @param configfile The path and name of the dynamic configuration file to read in. @param data A reference to the transient data class used to pass information between Tools. + bool Execute(); ///< Execute function used to perform Tool purpose. + bool Finalise(); ///< Finalise funciton used to clean up resources. + + bool does_file_exist(const char *fileName); ///< Function to check if a certain file exists + void draw_red_box(TH2 *hist, int bin_x, int bin_y, int verbosity); ///< Function to draw rectangles around PMTs that seem to be misaligned + + private: + + //define input variables + std::string outputfile; + std::string stabilityfile; + std::string geometryfile; + double diffuser_x; + double diffuser_y; + double diffuser_z; + double tolerance_charge; + double tolerance_time; + std::string FitMethod; + bool use_tapplication; + int verbose; + + //define ANNIEEvent variables + int evnum; + int runnumber; + TimeClass* EventTime=nullptr; + std::vector* TriggerData; + BeamStatusClass* BeamStatus=nullptr; + std::map>* MCHits = nullptr; + Geometry *geom = nullptr; + std::vector pmt_detkeys; + + + //define useful variables + const double n_water = 1.33; + const double c_vacuum= 2.99792E8; //in m/s + int n_tank_pmts, n_mrd_pmts, n_veto_pmts, n_lappds; + double tank_center_x, tank_center_y, tank_center_z; + std::map map_type_radius; //stores the radius information of the respective PMT as a function of the PMT type + + //define monitoring variables + std::map bad_time, bad_charge; + + //define arrays for storing PMT calibration values + std::map x_PMT; + std::map y_PMT; + std::map z_PMT; + std::map rho_PMT; + std::map phi_PMT; + std::map radius_PMT; + std::map PMT_ishit; + std::map mean_charge_fit; + std::map mean_time_fit; + std::map rms_charge_fit; + std::map rms_time_fit; + std::map expected_time; + + //define histograms for stability plots + TH1F *hist_tubeid = nullptr; + TH1F *hist_charge = nullptr; + TH1F *hist_time = nullptr; + std::map hist_charge_singletube; + std::map hist_time_singletube; + TH1F *hist_charge_mean = nullptr; + TH1F *hist_time_mean = nullptr; + TH1F *hist_time_dev = nullptr; + TH2F *hist_charge_2D_y_phi = nullptr; + TH2F *hist_time_2D_y_phi = nullptr; + TH2F *hist_time_2D_y_phi_mean = nullptr; + + //container for red boxes surrounding badly calibrated PMTs + std::vector vector_tbox; + + //container for TF1s used to fit the distributions + std::vector vector_tf1; + + //define graphs for stability plots + TGraphErrors *gr_stability = nullptr; + TGraphErrors *gr_stability_time = nullptr; + TGraphErrors *gr_stability_time_mean = nullptr; + + //define file to save data + TFile *file_out = nullptr; + TFile *help_file = nullptr; + + //define overview canvas + TCanvas *canvas_overview = nullptr; + TCanvas *canvas_overview2 = nullptr; + + //define TApplication for possibility of interactive plotting + TApplication *app_stability = nullptr; }; diff --git a/configfiles/TankCalibrationDiffuser/DummyToolConfig b/configfiles/TankCalibrationDiffuser/DummyToolConfig deleted file mode 100644 index 95cad88d2..000000000 --- a/configfiles/TankCalibrationDiffuser/DummyToolConfig +++ /dev/null @@ -1,3 +0,0 @@ -# Dummy config file - -verbose 2 \ No newline at end of file diff --git a/configfiles/TankCalibrationDiffuser/LoadWCSimConfig b/configfiles/TankCalibrationDiffuser/LoadWCSimConfig index 22f77bce4..b24acd8ab 100644 --- a/configfiles/TankCalibrationDiffuser/LoadWCSimConfig +++ b/configfiles/TankCalibrationDiffuser/LoadWCSimConfig @@ -1,6 +1,7 @@ # LoadWCSim Config File InputFile /ANNIECode/simulations/LaserCalibrationFiles/FakeRun_1000photons_5000events_timeshift_10.root +#InputFile /ANNIECode/simulations/LaserCalibrationFiles/wcsim_photonbomb_1000_5000events_ANNIEp2v6.root WCSimVersion 3 ## should reflect the WCSim version of the files being loaded HistoricTriggeroffset 0 ## time offset of digits relative to the trigger UseDigitSmearedTime 1 ## whether to use smeared digit time (T), or true time of first photon (F) diff --git a/configfiles/TankCalibrationDiffuser/TankCalibrationDiffuserConfig b/configfiles/TankCalibrationDiffuser/TankCalibrationDiffuserConfig index 9432c2228..5cc5f59be 100644 --- a/configfiles/TankCalibrationDiffuser/TankCalibrationDiffuserConfig +++ b/configfiles/TankCalibrationDiffuser/TankCalibrationDiffuserConfig @@ -1,13 +1,12 @@ # TankCalibrationDiffuser Config File -OutputFile /ANNIECode/ToolAnalysis/histogrammed_pmt_charge.root #Output root file for the current calibration run -StabilityFile /ANNIECode/ToolAnalysis/stability_pmtcharge.root #Output root file that displays the stability over multiple runs -GeometryFile configfiles/TankCalibrationDiffuser/geofile_128PMTs.txt #Geometry file specifying the different radii of the PMTs (taken from file geofile.txt in WCSim installation directory for the respective installation) -DiffuserX 0. #x-position of the diffuser ball -DiffuserY 0. #y-position of the diffuser ball -DiffuserZ 0. #z-position of the diffuser ball +OutputFile SimulationCalibrationTest2 #Output root prefix name for the current run +DiffuserX 0. #x-position of the diffuser ball +DiffuserY 0. #y-position of the diffuser ball +DiffuserZ 0. #z-position of the diffuser ball ToleranceCharge 0.5 #tolerance of fit single p.e. value for being classified as a bad PMT ToleranceTime 0.5 #tolerance of mean time value [ns] for being classified as a bad PMT +FitMethod Gaus2Exp #fit function for charge, options: Gaus2Exp (2 times gaus + exp), Gaus2 (2 times gaus), Gaus (single gaus) TApplication 0 #0/1, depending on whether plots should be shown interactively or not verbose 1 #verbosity of the application diff --git a/configfiles/TankCalibrationDiffuser/geofile_128PMTs.txt b/configfiles/TankCalibrationDiffuser/geofile_128PMTs.txt deleted file mode 100644 index 10034a200..000000000 --- a/configfiles/TankCalibrationDiffuser/geofile_128PMTs.txt +++ /dev/null @@ -1,128 +0,0 @@ -R5912HQE -R7081 -R7081 -R7081 -R5912HQE -R5912HQE -R7081 -R5912HQE -R5912HQE -R7081 -R7081 -R7081 -R5912HQE -R5912HQE -R7081 -R5912HQE -R5912HQE -R7081 -R7081 -R7081 -R5912HQE -R5912HQE -R7081 -R5912HQE -R5912HQE -R7081 -R7081 -R7081 -R5912HQE -R5912HQE -R7081 -R5912HQE -R5912HQE -R7081 -R7081 -R7081 -R5912HQE -R5912HQE -R7081 -R5912HQE -R5912HQE -R7081 -R7081 -R7081 -R5912HQE -R5912HQE -R7081 -R5912HQE -R5912HQE -R7081 -R7081 -R7081 -R5912HQE -R5912HQE -R7081 -R5912HQE -R5912HQE -R7081 -R7081 -R7081 -R5912HQE -R5912HQE -R7081 -R5912HQE -R5912HQE -R5912HQE -R5912HQE -R5912HQE -R5912HQE -R5912HQE -R5912HQE -R5912HQE -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB diff --git a/configfiles/TankCalibrationDiffuser/geofile_143PMTs.txt b/configfiles/TankCalibrationDiffuser/geofile_143PMTs.txt deleted file mode 100644 index 52f122c83..000000000 --- a/configfiles/TankCalibrationDiffuser/geofile_143PMTs.txt +++ /dev/null @@ -1,143 +0,0 @@ -R5912HQE -R7081 -R7081 -R7081HQE -R7081 -R5912HQE -R7081 -R5912HQE -R7081 -R7081 -R5912HQE -R5912HQE -R7081 -R7081 -R7081HQE -R7081 -R5912HQE -R5912HQE -R7081 -R5912HQE -R7081 -R7081 -R5912HQE -R5912HQE -R5912HQE -R7081 -R7081HQE -R7081 -R5912HQE -R5912HQE -R7081 -R5912HQE -R7081 -R7081 -R5912HQE -R5912HQE -R5912HQE -R7081 -R7081HQE -R7081 -R5912HQE -R5912HQE -R7081 -R7081HQE -R7081 -R7081 -R5912HQE -R5912HQE -R5912HQE -R7081 -R7081HQE -R7081 -R5912HQE -R5912HQE -R7081 -R7081HQE -R7081 -R7081 -R5912HQE -R5912HQE -R7081 -R7081 -R7081HQE -R7081 -R5912HQE -R5912HQE -R7081 -R5912HQE -R7081 -R7081 -R5912HQE -R5912HQE -R7081 -R7081 -R7081HQE -R7081 -R5912HQE -R5912HQE -R7081 -R5912HQE -R7081 -R7081 -R5912HQE -R5912HQE -R7081 -R7081 -R7081HQE -R7081 -R5912HQE -R5912HQE -R7081 -R5912HQE -R7081 -R7081 -R5912HQE -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -D784KFLB -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -R7081 -EMI9954KB -EMI9954KB -EMI9954KB -EMI9954KB -EMI9954KB -EMI9954KB -EMI9954KB -EMI9954KB \ No newline at end of file