diff --git a/src/plugins/Analysis/SConscript b/src/plugins/Analysis/SConscript index bfc1877d10..e1c4d1fb9c 100644 --- a/src/plugins/Analysis/SConscript +++ b/src/plugins/Analysis/SConscript @@ -7,23 +7,25 @@ Import('*') subdirs = ['ReactionFilter', 'b1pi_hists', 'ReactionEfficiency'] subdirs.extend(['DAQTree', 'DAQTreeBCAL']) subdirs.extend(['monitoring_hists']) -subdirs.extend(['p2pi_hists','p3pi_hists','p4pi_hists','p2k_hists']) -subdirs.extend(['p2gamma_hists','ppi0gamma_hists','p2pi0_hists']) -subdirs.extend(['npp_hists','cpp_hists']) +subdirs.extend(['p2pi_hists','p3pi_hists','ppi0gamma_hists']) subdirs.extend(['fcal_charged','TPOL_tree','mcthrown_tree']) -subdirs.extend(['trackeff_hists','imaging']) -subdirs.extend(['pid_dirc','truth_dirc', 'lut_dirc', 'dirc_hists', 'dirc_reactions', 'dirc_tree', 'TRD_hists']) -subdirs.extend(['B3pi_eff_missgamma','B3pi_eff_misspip','B3pi_eff_misspim','B3pi_eff_missprot', 'mcthrown_hists']) -subdirs.extend(['compton']) -subdirs.extend(['src-ct']) +subdirs.extend(['trackeff_hists','imaging','dirc_reactions','TRD_hists']) SConscript(dirs=subdirs, exports='env osname', duplicate=0) # Optional targets -optdirs = ['phys_tree', 'pedestals','bcal_calib','bcal_calib_cosmic_cdc', 'CPPMVAtree'] - -optdirs.extend(['acceptance_hists', 'mcthrown_hists', 'cdc_hists', 'fcal_hists']) -optdirs.extend(['dc_alignment','p2pi_trees','Z2pi_trees','ccal_hits','ccal_display','fcal_led','F250_mode10_pedestal','Z2pi0_trees']) - +optdirs = ['pedestals','bcal_calib','bcal_calib_cosmic_cdc', 'CPPMVAtree'] + +optdirs.extend(['mcthrown_hists', 'cdc_hists', 'fcal_hists']) +optdirs.extend(['p4pi_hists','p2k_hists','p2gamma_hists','p2pi0_hists']) +optdirs.extend(['B3pi_eff_missgamma','B3pi_eff_misspip','B3pi_eff_misspim','B3pi_eff_missprot']) +subdirs.extend(['pid_dirc','truth_dirc', 'lut_dirc', 'dirc_hists', 'dirc_tree']) +optdirs.extend(['dc_alignment','p2pi_trees','fcal_led','F250_mode10_pedestal']) +optdirs.extend(['npp_hists','cpp_hists']) +optdirs.extend(['Z2pi0_trees','Z2pi_trees','ccal_hits','ccal_display']) +optdirs.extend(['compton','src-ct']) + +optdirs.extend(['retired']) + sbms.OptionallyBuild(env, optdirs) diff --git a/src/plugins/Analysis/acceptance_hists/DEventProcessor_acceptance_hists.cc b/src/plugins/Analysis/retired/acceptance_hists/DEventProcessor_acceptance_hists.cc similarity index 100% rename from src/plugins/Analysis/acceptance_hists/DEventProcessor_acceptance_hists.cc rename to src/plugins/Analysis/retired/acceptance_hists/DEventProcessor_acceptance_hists.cc diff --git a/src/plugins/Analysis/acceptance_hists/DEventProcessor_acceptance_hists.h b/src/plugins/Analysis/retired/acceptance_hists/DEventProcessor_acceptance_hists.h similarity index 100% rename from src/plugins/Analysis/acceptance_hists/DEventProcessor_acceptance_hists.h rename to src/plugins/Analysis/retired/acceptance_hists/DEventProcessor_acceptance_hists.h diff --git a/src/plugins/Analysis/acceptance_hists/Makefile b/src/plugins/Analysis/retired/acceptance_hists/Makefile similarity index 100% rename from src/plugins/Analysis/acceptance_hists/Makefile rename to src/plugins/Analysis/retired/acceptance_hists/Makefile diff --git a/src/plugins/Analysis/acceptance_hists/SConscript b/src/plugins/Analysis/retired/acceptance_hists/SConscript similarity index 100% rename from src/plugins/Analysis/acceptance_hists/SConscript rename to src/plugins/Analysis/retired/acceptance_hists/SConscript diff --git a/src/plugins/Analysis/phys_tree/DEventProcessor_phys_tree.cc b/src/plugins/Analysis/retired/phys_tree/DEventProcessor_phys_tree.cc similarity index 100% rename from src/plugins/Analysis/phys_tree/DEventProcessor_phys_tree.cc rename to src/plugins/Analysis/retired/phys_tree/DEventProcessor_phys_tree.cc diff --git a/src/plugins/Analysis/phys_tree/DEventProcessor_phys_tree.h b/src/plugins/Analysis/retired/phys_tree/DEventProcessor_phys_tree.h similarity index 100% rename from src/plugins/Analysis/phys_tree/DEventProcessor_phys_tree.h rename to src/plugins/Analysis/retired/phys_tree/DEventProcessor_phys_tree.h diff --git a/src/plugins/Analysis/phys_tree/Event.cc b/src/plugins/Analysis/retired/phys_tree/Event.cc similarity index 100% rename from src/plugins/Analysis/phys_tree/Event.cc rename to src/plugins/Analysis/retired/phys_tree/Event.cc diff --git a/src/plugins/Analysis/phys_tree/Event.h b/src/plugins/Analysis/retired/phys_tree/Event.h similarity index 100% rename from src/plugins/Analysis/phys_tree/Event.h rename to src/plugins/Analysis/retired/phys_tree/Event.h diff --git a/src/plugins/Analysis/phys_tree/Makefile b/src/plugins/Analysis/retired/phys_tree/Makefile similarity index 100% rename from src/plugins/Analysis/phys_tree/Makefile rename to src/plugins/Analysis/retired/phys_tree/Makefile diff --git a/src/plugins/Analysis/phys_tree/Particle.h b/src/plugins/Analysis/retired/phys_tree/Particle.h similarity index 100% rename from src/plugins/Analysis/phys_tree/Particle.h rename to src/plugins/Analysis/retired/phys_tree/Particle.h diff --git a/src/plugins/Analysis/phys_tree/SConscript b/src/plugins/Analysis/retired/phys_tree/SConscript similarity index 100% rename from src/plugins/Analysis/phys_tree/SConscript rename to src/plugins/Analysis/retired/phys_tree/SConscript diff --git a/src/plugins/Calibration/SConscript b/src/plugins/Calibration/SConscript index 199c1f2881..23890a5026 100644 --- a/src/plugins/Calibration/SConscript +++ b/src/plugins/Calibration/SConscript @@ -5,8 +5,11 @@ Import('*') # update this when plugins are added -subdirs = ['BCAL_SiPM_saturation','BCAL_attenlength_gainratio','BCAL_gainmatrix','BCAL_point_time','BCAL_point_calib','BCAL_saturation','BCAL_TDC_Timing','CDC_amp','HLDetectorTiming','TAGH_timewalk','CDC_TimeToDistance','PS_timing','PSC_TW','PS_E_calib','st_tw_corr_auto','ST_Propagation_Time','ST_Tresolution','TAGM_TW','TOF_calib','FCALLEDTree','FCAL_LED_shifts','FCAL_TimingOffsets_Primex'] +subdirs = ['BCAL_attenlength_gainratio','BCAL_gainmatrix','BCAL_TDC_Timing','CDC_amp','HLDetectorTiming','TAGH_timewalk','CDC_TimeToDistance','PS_timing','PSC_TW','PS_E_calib','st_tw_corr_auto','ST_Propagation_Time','ST_Tresolution','TAGM_TW','TOF_calib','FCALLEDTree','FCAL_TimingOffsets_Primex'] SConscript(dirs=subdirs, exports='env osname', duplicate=0) -sbms.OptionallyBuild(env, ['FCALgains','FCALpedestals','FCALpulsepeak','FCAL_TimingOffsets','FCAL_Pi0HFA','FCAL_Pi0TOF']) +optdirs = ['FCALpedestals','FCALpulsepeak','FCAL_TimingOffsets','FCAL_Pi0TOF'] +optdirs.extend(['BCAL_saturation','BCAL_SiPM_saturation','FCAL_LED_shifts']) + +sbms.OptionallyBuild(env, optdirs) diff --git a/src/plugins/Calibration/BCAL_point_calib/JEventProcessor_BCAL_point_calib.cc b/src/plugins/Calibration/retired/BCAL_point_calib/JEventProcessor_BCAL_point_calib.cc similarity index 100% rename from src/plugins/Calibration/BCAL_point_calib/JEventProcessor_BCAL_point_calib.cc rename to src/plugins/Calibration/retired/BCAL_point_calib/JEventProcessor_BCAL_point_calib.cc diff --git a/src/plugins/Calibration/BCAL_point_calib/JEventProcessor_BCAL_point_calib.h b/src/plugins/Calibration/retired/BCAL_point_calib/JEventProcessor_BCAL_point_calib.h similarity index 100% rename from src/plugins/Calibration/BCAL_point_calib/JEventProcessor_BCAL_point_calib.h rename to src/plugins/Calibration/retired/BCAL_point_calib/JEventProcessor_BCAL_point_calib.h diff --git a/src/plugins/Calibration/BCAL_point_calib/SConscript b/src/plugins/Calibration/retired/BCAL_point_calib/SConscript similarity index 100% rename from src/plugins/Calibration/BCAL_point_calib/SConscript rename to src/plugins/Calibration/retired/BCAL_point_calib/SConscript diff --git a/src/plugins/Calibration/BCAL_point_calib/c_eff_linear_fit.C b/src/plugins/Calibration/retired/BCAL_point_calib/c_eff_linear_fit.C similarity index 100% rename from src/plugins/Calibration/BCAL_point_calib/c_eff_linear_fit.C rename to src/plugins/Calibration/retired/BCAL_point_calib/c_eff_linear_fit.C diff --git a/src/plugins/Calibration/BCAL_point_calib/c_eff_pol2.C b/src/plugins/Calibration/retired/BCAL_point_calib/c_eff_pol2.C similarity index 100% rename from src/plugins/Calibration/BCAL_point_calib/c_eff_pol2.C rename to src/plugins/Calibration/retired/BCAL_point_calib/c_eff_pol2.C diff --git a/src/plugins/Calibration/BCAL_point_calib/z_point_pol2.C b/src/plugins/Calibration/retired/BCAL_point_calib/z_point_pol2.C similarity index 100% rename from src/plugins/Calibration/BCAL_point_calib/z_point_pol2.C rename to src/plugins/Calibration/retired/BCAL_point_calib/z_point_pol2.C diff --git a/src/plugins/Calibration/BCAL_point_time/JEventProcessor_BCAL_point_time.cc b/src/plugins/Calibration/retired/BCAL_point_time/JEventProcessor_BCAL_point_time.cc similarity index 100% rename from src/plugins/Calibration/BCAL_point_time/JEventProcessor_BCAL_point_time.cc rename to src/plugins/Calibration/retired/BCAL_point_time/JEventProcessor_BCAL_point_time.cc diff --git a/src/plugins/Calibration/BCAL_point_time/JEventProcessor_BCAL_point_time.h b/src/plugins/Calibration/retired/BCAL_point_time/JEventProcessor_BCAL_point_time.h similarity index 100% rename from src/plugins/Calibration/BCAL_point_time/JEventProcessor_BCAL_point_time.h rename to src/plugins/Calibration/retired/BCAL_point_time/JEventProcessor_BCAL_point_time.h diff --git a/src/plugins/Calibration/BCAL_point_time/SConscript b/src/plugins/Calibration/retired/BCAL_point_time/SConscript similarity index 100% rename from src/plugins/Calibration/BCAL_point_time/SConscript rename to src/plugins/Calibration/retired/BCAL_point_time/SConscript diff --git a/src/plugins/Calibration/BCAL_point_time/macros/plot_BCAL_point_time.C b/src/plugins/Calibration/retired/BCAL_point_time/macros/plot_BCAL_point_time.C similarity index 100% rename from src/plugins/Calibration/BCAL_point_time/macros/plot_BCAL_point_time.C rename to src/plugins/Calibration/retired/BCAL_point_time/macros/plot_BCAL_point_time.C diff --git a/src/plugins/Calibration/FCAL_Pi0HFA/FitScripts/FitGains.C b/src/plugins/Calibration/retired/FCAL_Pi0HFA/FitScripts/FitGains.C similarity index 100% rename from src/plugins/Calibration/FCAL_Pi0HFA/FitScripts/FitGains.C rename to src/plugins/Calibration/retired/FCAL_Pi0HFA/FitScripts/FitGains.C diff --git a/src/plugins/Calibration/FCAL_Pi0HFA/JEventProcessor_FCAL_Pi0HFA.cc b/src/plugins/Calibration/retired/FCAL_Pi0HFA/JEventProcessor_FCAL_Pi0HFA.cc similarity index 100% rename from src/plugins/Calibration/FCAL_Pi0HFA/JEventProcessor_FCAL_Pi0HFA.cc rename to src/plugins/Calibration/retired/FCAL_Pi0HFA/JEventProcessor_FCAL_Pi0HFA.cc diff --git a/src/plugins/Calibration/FCAL_Pi0HFA/JEventProcessor_FCAL_Pi0HFA.h b/src/plugins/Calibration/retired/FCAL_Pi0HFA/JEventProcessor_FCAL_Pi0HFA.h similarity index 100% rename from src/plugins/Calibration/FCAL_Pi0HFA/JEventProcessor_FCAL_Pi0HFA.h rename to src/plugins/Calibration/retired/FCAL_Pi0HFA/JEventProcessor_FCAL_Pi0HFA.h diff --git a/src/plugins/Calibration/FCAL_Pi0HFA/SConscript b/src/plugins/Calibration/retired/FCAL_Pi0HFA/SConscript similarity index 100% rename from src/plugins/Calibration/FCAL_Pi0HFA/SConscript rename to src/plugins/Calibration/retired/FCAL_Pi0HFA/SConscript diff --git a/src/plugins/Calibration/FCALgains/JEventProcessor_FCALgains.cc b/src/plugins/Calibration/retired/FCALgains/JEventProcessor_FCALgains.cc similarity index 100% rename from src/plugins/Calibration/FCALgains/JEventProcessor_FCALgains.cc rename to src/plugins/Calibration/retired/FCALgains/JEventProcessor_FCALgains.cc diff --git a/src/plugins/Calibration/FCALgains/JEventProcessor_FCALgains.h b/src/plugins/Calibration/retired/FCALgains/JEventProcessor_FCALgains.h similarity index 100% rename from src/plugins/Calibration/FCALgains/JEventProcessor_FCALgains.h rename to src/plugins/Calibration/retired/FCALgains/JEventProcessor_FCALgains.h diff --git a/src/plugins/Calibration/FCALgains/SConscript b/src/plugins/Calibration/retired/FCALgains/SConscript similarity index 100% rename from src/plugins/Calibration/FCALgains/SConscript rename to src/plugins/Calibration/retired/FCALgains/SConscript diff --git a/src/plugins/SConscript b/src/plugins/SConscript index f79730c83e..c9057a6174 100644 --- a/src/plugins/SConscript +++ b/src/plugins/SConscript @@ -5,6 +5,8 @@ Import('*') Import('env osname') -subdirs = ['Analysis', 'Utilities', 'monitoring', 'Calibration', 'include', 'Alignment'] +subdirs = ['Analysis', 'Utilities', 'monitoring', 'Calibration', 'include'] +sbms.OptionallyBuild(env, ['Alignment']) + SConscript(dirs=subdirs, exports='env osname', duplicate=0) diff --git a/src/plugins/Utilities/SConscript b/src/plugins/Utilities/SConscript index 3a10df794e..dea1c9438a 100644 --- a/src/plugins/Utilities/SConscript +++ b/src/plugins/Utilities/SConscript @@ -5,14 +5,18 @@ Import('*') # Default targets (always built) subdirs = ['danarest', 'evio_writer', 'evio-hddm'] -subdirs.extend(['2trackskim', 'pi0bcalskim', 'pi0fcalskim', 'twogamma_fcal_skim', 'run_summary', 'track_skimmer', 'trackeff_missing','ps_skim', 'trigger_skims', 'bigevents_skim', 'coherent_peak_skim','exclusivepi0skim','randomtrigger_skim','pi0fcaltofskim','single_neutral_skim','compton_neutral_skim','eta2g_primexd_skim','eta6g_primexd_skim','etapi0_primexd_skim', 'cdcbcal_skim', 'cdc_goodtrack_skim']) - -subdirs.extend(['Pi0Finder', 'EventTagPi0','es_test','omega_skim','cal_high_energy_skim', 'syncskim', 'dedx_tree', 'phi_skim', 'npp_skim', 'cpp_skim']) +subdirs.extend(['pi0bcalskim', 'trackeff_missing','ps_skim', 'trigger_skims', 'exclusivepi0skim', 'randomtrigger_skim', 'pi0fcaltofskim' ]) +subdirs.extend(['omega_skim', 'syncskim']) SConscript(dirs=subdirs, exports='env osname', duplicate=0) # Optional targets -optdirs = ['danahddm', 'dumpcandidates', 'dumpthrowns', 'l3bdt'] -optdirs.extend(['merge_rawevents', 'syncskim', 'DAQ', 'TTab', 'rawevent']) +optdirs = [ '2trackskim', 'bigevents_skim', 'coherent_peak_skim' , 'twogamma_fcal_skim' ] +optdirs.extend(['pi0fcalskim', 'single_neutral_skim','compton_neutral_skim' ]) +optdirs.extend(['eta2g_primexd_skim','eta6g_primexd_skim','etapi0_primexd_skim' ]) +optdirs.extend(['cdcbcal_skim', 'cdc_goodtrack_skim']) +optdirs.extend(['run_summary', 'track_skimmer' ]) +optdirs.extend(['merge_rawevents', 'syncskim']) +optdirs.extend(['Pi0Finder', 'EventTagPi0','cal_high_energy_skim', 'phi_skim', 'npp_skim', 'cpp_skim']) optdirs.extend(['cdc_amp_t','cdc_echo','cdc_emu','cdc_scan', 'fmwpc_scan', 'epem_ml_skim']) sbms.OptionallyBuild(env, optdirs) diff --git a/src/plugins/Utilities/danahddm/JEventProcessor_danahddm.cc b/src/plugins/Utilities/retired/danahddm/JEventProcessor_danahddm.cc similarity index 100% rename from src/plugins/Utilities/danahddm/JEventProcessor_danahddm.cc rename to src/plugins/Utilities/retired/danahddm/JEventProcessor_danahddm.cc diff --git a/src/plugins/Utilities/danahddm/JEventProcessor_danahddm.h b/src/plugins/Utilities/retired/danahddm/JEventProcessor_danahddm.h similarity index 100% rename from src/plugins/Utilities/danahddm/JEventProcessor_danahddm.h rename to src/plugins/Utilities/retired/danahddm/JEventProcessor_danahddm.h diff --git a/src/plugins/Utilities/danahddm/Makefile b/src/plugins/Utilities/retired/danahddm/Makefile similarity index 100% rename from src/plugins/Utilities/danahddm/Makefile rename to src/plugins/Utilities/retired/danahddm/Makefile diff --git a/src/plugins/Utilities/danahddm/SConscript b/src/plugins/Utilities/retired/danahddm/SConscript similarity index 100% rename from src/plugins/Utilities/danahddm/SConscript rename to src/plugins/Utilities/retired/danahddm/SConscript diff --git a/src/plugins/Utilities/retired/dedx_tree/JEventProcessor_dedx_tree.cc b/src/plugins/Utilities/retired/dedx_tree/JEventProcessor_dedx_tree.cc new file mode 100644 index 0000000000..917badb7b8 --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/JEventProcessor_dedx_tree.cc @@ -0,0 +1,199 @@ +// $Id$ +// +// File: JEventProcessor_dedx_tree.cc +// Created: Wed Aug 28 14:48:25 EDT 2019 +// Creator: njarvis (on Linux albert.phys.cmu.edu 3.10.0-693.5.2.el7.x86_64 x86_64) +// + +#include "JEventProcessor_dedx_tree.h" + + +// Routine used to create our JEventProcessor +#include +#include + +extern "C"{ +void InitPlugin(JApplication *app){ + InitJANAPlugin(app); + app->Add(new JEventProcessor_dedx_tree()); +} +} // "C" + +//define static local variable //declared in header file +thread_local DTreeFillData JEventProcessor_dedx_tree::dTreeFillData; + +//------------------ +// JEventProcessor_dedx_tree (Constructor) +//------------------ +JEventProcessor_dedx_tree::JEventProcessor_dedx_tree() +{ + SetTypeName("JEventProcessor_dedx_tree"); +} + +//------------------ +// ~JEventProcessor_dedx_tree (Destructor) +//------------------ +JEventProcessor_dedx_tree::~JEventProcessor_dedx_tree() +{ + +} + +//------------------ +// init +//------------------ +void JEventProcessor_dedx_tree::Init() +{ + // This is called once at program startup. + + //TTREE INTERFACE + //MUST DELETE WHEN FINISHED: OR ELSE DATA WON'T BE SAVED!!! + dTreeInterface = DTreeInterface::Create_DTreeInterface("dedx_tree", "dedx_tree.root"); + + //TTREE BRANCHES + DTreeBranchRegister locTreeBranchRegister; + + locTreeBranchRegister.Register_Single("eventnumber"); + locTreeBranchRegister.Register_Single("ctracks"); + + locTreeBranchRegister.Register_Single("ndedxhits"); + + locTreeBranchRegister.Register_Single("x"); + locTreeBranchRegister.Register_Single("y"); + locTreeBranchRegister.Register_Single("z"); + locTreeBranchRegister.Register_Single("r"); + + locTreeBranchRegister.Register_Single("phi"); + locTreeBranchRegister.Register_Single("theta"); + + locTreeBranchRegister.Register_Single("charge"); + locTreeBranchRegister.Register_Single("p"); + locTreeBranchRegister.Register_Single("dedx"); + locTreeBranchRegister.Register_Single("dedx_int"); + locTreeBranchRegister.Register_Single("dedx_corr"); + locTreeBranchRegister.Register_Single("dedx_int_corr"); + + //REGISTER BRANCHES + dTreeInterface->Create_Branches(locTreeBranchRegister); +} + +//------------------ +// brun +//------------------ +void JEventProcessor_dedx_tree::BeginRun(const std::shared_ptr &event) +{ + // This is called whenever the run number changes +} + +//------------------ +// evnt +//------------------ +void JEventProcessor_dedx_tree::Process(const std::shared_ptr &event) +{ + // This is called for every event. Use of common resources like writing + // to a file or filling a histogram should be mutex protected. Using + // event->Get(...) to get reconstructed objects (and thereby activating the + // reconstruction algorithm) should be done outside of any mutex lock + // since multiple threads may call this method at the same time. + // Here's an example: + // + // vector mydataclasses; + // event->Get(mydataclasses); + // + // japp->RootFillLock(this); + // ... fill historgrams or trees ... + // japp->RootFillUnLock(this); + + auto eventnumber = event->GetEventNumber(); + + // select events with physics events, i.e., not LED and other front panel triggers + const DTrigger* locTrigger = NULL; + event->GetSingle(locTrigger); + if(locTrigger->Get_L1FrontPanelTriggerBits() != 0) return; + + const DVertex* locVertex = NULL; + event->GetSingle(locVertex); + double z = locVertex->dSpacetimeVertex.Z(); + if ((z < 45.0) || (z > 85.0)) return; + + double x = locVertex->dSpacetimeVertex.X(); + double y = locVertex->dSpacetimeVertex.Y(); + double r = sqrt (x*x + y*y); + + + vector ctracks; + event->Get(ctracks); + + if ((int)ctracks.size() ==0) return; + + dTreeFillData.Fill_Single("eventnumber",(ULong64_t)eventnumber); + + dTreeFillData.Fill_Single("ctracks",(UInt_t)ctracks.size()); + + + for (uint32_t i=0; iGet_BestFOM(); + + if (hyp != NULL) { + + const DTrackTimeBased *track = hyp->Get_TrackTimeBased(); + + int nhits = (int)track->dNumHitsUsedFordEdx_CDC; + + if (nhits==0) continue; + + dTreeFillData.Fill_Single("ndedxhits",(UInt_t)nhits); + + double charge = track->charge(); + DVector3 mom = track->momentum(); + double p = mom.Mag(); + + double theta_degrees = mom.Theta() * 180.0/3.14159; + double phi_degrees = mom.Phi() * 180.0/3.14159; + + double dedx = 1.0e6*track->ddEdx_CDC_amp; + double dedxfromintegral = 1.0e6*track->ddEdx_CDC; + double dedx_corr = 1.0e6*hyp->Get_dEdx_CDC_amp(); + double dedxfromintegral_corr = 1.0e6*hyp->Get_dEdx_CDC_int(); + + dTreeFillData.Fill_Single("x",x); + dTreeFillData.Fill_Single("y",y); + dTreeFillData.Fill_Single("z",z); + dTreeFillData.Fill_Single("r",r); + dTreeFillData.Fill_Single("phi",phi_degrees); + dTreeFillData.Fill_Single("theta",theta_degrees); + + dTreeFillData.Fill_Single("charge",(Int_t)charge); + dTreeFillData.Fill_Single("p",p); + dTreeFillData.Fill_Single("dedx",dedx); + dTreeFillData.Fill_Single("dedx_int",dedxfromintegral); + dTreeFillData.Fill_Single("dedx_corr",dedx_corr); + dTreeFillData.Fill_Single("dedx_int_corr",dedxfromintegral_corr); + + //FILL TTREE + dTreeInterface->Fill(dTreeFillData); + } + + } +} + +//------------------ +// erun +//------------------ +void JEventProcessor_dedx_tree::EndRun() +{ + // This is called whenever the run number changes, before it is + // changed to give you a chance to clean up before processing + // events from the next run number. +} + +//------------------ +// fini +//------------------ +void JEventProcessor_dedx_tree::Finish() +{ + // Called before program exit after event processing is finished. + + delete dTreeInterface; //saves trees to file, closes file +} + diff --git a/src/plugins/Utilities/retired/dedx_tree/JEventProcessor_dedx_tree.h b/src/plugins/Utilities/retired/dedx_tree/JEventProcessor_dedx_tree.h new file mode 100644 index 0000000000..beb27ea24d --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/JEventProcessor_dedx_tree.h @@ -0,0 +1,49 @@ +// $Id$ +// +// File: JEventProcessor_dedx_tree.h +// Created: Wed Aug 28 14:48:25 EDT 2019 +// Creator: njarvis (on Linux albert.phys.cmu.edu 3.10.0-693.5.2.el7.x86_64 x86_64) +// + +#ifndef _JEventProcessor_dedx_tree_ +#define _JEventProcessor_dedx_tree_ + +#include +#include + +#include "TRIGGER/DTrigger.h" +#include "TRACKING/DTrackTimeBased.h" +#include "PID/DVertex.h" +#include "PID/DChargedTrack.h" +#include "PID/DChargedTrackHypothesis.h" + +#include "ANALYSIS/DTreeInterface.h" + +#include + + +class JEventProcessor_dedx_tree:public JEventProcessor{ + public: + JEventProcessor_dedx_tree(); + ~JEventProcessor_dedx_tree(); + + private: + void Init() override; + void BeginRun(const std::shared_ptr& event) override; + void Process(const std::shared_ptr& event) override; + void EndRun() override; + void Finish() override; + + //TREE + DTreeInterface* dTreeInterface; + + //thread_local: Each thread has its own object: no lock needed + //important: manages it's own data internally: don't want to call new/delete every event! + + static thread_local DTreeFillData dTreeFillData; + + +}; + +#endif // _JEventProcessor_dedx_tree_ + diff --git a/src/plugins/Utilities/retired/dedx_tree/README.md b/src/plugins/Utilities/retired/dedx_tree/README.md new file mode 100644 index 0000000000..e13e6c51f3 --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/README.md @@ -0,0 +1,53 @@ +# CDC dE/dx correction procedure for variation with theta + +The first iteration of this process is documented in [GlueX-doc-4693](https://halldweb.jlab.org/doc-private/DocDB/ShowDocument?docid=4693). +The scripts are stored in the subdirectory scripts. +The correction matrix is loaded into the CCDB table /CDC/dedx\_theta. + +## Prerequisites: +* 300+ evio files processed with the dedx\_tree plugin, with dE/dx truncation and the dE/dx correction switched off ( -PPID:CDC_CORRECT_DEDX_THETA=0 and -PPID:CDC_TRUNCATE_DEDX=0 ) + +* Simulated data processed similarly. + + +## Process overview: + +1. makebinnedhistos.C - create many histograms of dedx vs p, for 1 degree bins in theta + +2. fitbinnedhistoslg.C - find the peak position for the proton and pion bands in each (p, theta) bin, using theta-specific functions to define the fit regions + +3. makesimsinglebinnedhistos.C - as (1), with simulated data + +4. fitbinnedhistoslg.C - as (2) with simulated data + +5. makescalefactors.C - calculate the scale factors to match the real data to the simulated data + +6. extendscalefactors.C - extrapolate the scale factors to populate the edges of a rectangular space in (dedx, theta) (root -b -q dedx_correction_factors.root extendscalefactors.C) + +7. smooth_histo.C - run a smoothing algorithm over the scale factors and write the matrix of scale factors into a text file which can be used in the reconstruction code. (root extended_dedx_correction_factors.root smooth_histo.C) + +8. checkmatrix.C - use the matrix contained in the text file to correct the data in the root trees and create a set of corrected binned histograms, which can be fitted as before to verify that the correction process produced acceptable results. + + +## Plotting scripts +* drawdedxcuts.C - plots functions on top of specified dEdx vs p histogram from histogram file +* drawalldedxcuts.C - loops through all histograms in histo file, saves 2x2 canvases of dedx functions over data +* plotdata.C - reads in fitresults.root or fittedsimresults.root, plots coloured bands of dedx peak position for one particle type +* plotf.C - reads in dedx_correction_factors.root, plots coloured bands of correction factors for one particle type +* plotfdedx.C - reads in dedx_correction_factors.root, plots dedx vs theta for the data used for the correction factors +* plot_cf_vs_dedx.C - reads in dedx_correction_factors.root, makes plots cf vs theta for 5 degree steps in theta +* dots.C - reads in dedx_correction_factors.root - makes a selection of plots of cf, dedx, theta including a 3D dot plot +* drawsurf2.C - use this with the tree file open, works with dedx_correction_factors.root or extended_dedx_correction_factors.root - draws surfaces, dots and surface+dots + + +## Other scripts +* make_histo_into_text_file.C - writes the histogram into a text file without smoothing + + + +A quick 3D plot of plot correction factors vs dedx, theta in colour scale can be made from the open root tree using +```sh +t->Draw("f:dedx:theta:f","","cont4 col") +``` + + diff --git a/src/plugins/Utilities/dumpcandidates/SConscript b/src/plugins/Utilities/retired/dedx_tree/SConscript similarity index 100% rename from src/plugins/Utilities/dumpcandidates/SConscript rename to src/plugins/Utilities/retired/dedx_tree/SConscript diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/checkmatrix.C b/src/plugins/Utilities/retired/dedx_tree/scripts/checkmatrix.C new file mode 100644 index 0000000000..b55cd9154b --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/checkmatrix.C @@ -0,0 +1,307 @@ +// this reads in a dedx_theta correction matrix +// reads in the dedx_tree root trees +// uses the matrix to correct dedx and writes out binned histograms into a new root file + +vector>CDC_DEDX_CORRECTION; +Float_t cdc_min_theta, cdc_max_theta; +Float_t cdc_min_dedx, cdc_max_dedx; +Float_t cdc_theta_step, cdc_dedx_step; +Int_t cdc_npoints_theta, cdc_npoints_dedx; + + +void Read_CDC_dEdx_CorrectionFile(void){ + + // FILE *dedxfile = fopen("/home/njarvis/phi/dsel/dedx_amp_corr_22July2020.txt","r"); + FILE *dedxfile = fopen("/home/njarvis/online/dedx/16hits/dedx_corrections_s5a_t_116_101.txt","r"); + fscanf(dedxfile,"%i values of theta\n",&cdc_npoints_theta); + fscanf(dedxfile,"%f min theta\n",&cdc_min_theta); + fscanf(dedxfile,"%f max theta\n",&cdc_max_theta); + fscanf(dedxfile,"%f theta step\n",&cdc_theta_step); + + fscanf(dedxfile,"%i values of dedx\n",&cdc_npoints_dedx); + + fscanf(dedxfile,"%f min dedx\n",&cdc_min_dedx); + fscanf(dedxfile,"%f max dedx\n",&cdc_max_dedx); + fscanf(dedxfile,"%f dedx step\n",&cdc_dedx_step); + fscanf(dedxfile,"\n"); + + vector dedx_cf_alltheta; + Float_t dedx_cf; + + // Store the scaling factors in vector>CDC_DEDX_CORRECTION; + + for (Int_t ii =0; ii= cdc_max_theta) { + thetabin1 = cdc_npoints_theta - 1; + thetabin2 = thetabin1; + } else { + thetabin1 = (Int_t)((theta_deg - cdc_min_theta)/cdc_theta_step); + thetabin2 = thetabin1 + 1; + } + + if (thisdedx <= cdc_min_dedx) { + dedxbin1 = 0; + dedxbin2 = dedxbin1; + + } else if (thisdedx >= cdc_max_dedx) { + dedxbin1 = cdc_npoints_dedx - 1; + dedxbin2 = dedxbin1; + } else { + dedxbin1 = (Int_t)((thisdedx - cdc_min_dedx)/cdc_dedx_step); + dedxbin2 = dedxbin1 + 1; + } + + Float_t dedxcf; + + if ((thetabin1 == thetabin2) && (dedxbin1 == dedxbin2)) { + + dedxcf = CDC_DEDX_CORRECTION[dedxbin1][thetabin1]; + + } else if (thetabin1 == thetabin2) { // interp dedx only + + Float_t cf1 = CDC_DEDX_CORRECTION[dedxbin1][thetabin1]; + Float_t cf2 = CDC_DEDX_CORRECTION[dedxbin2][thetabin1]; + + Float_t dedx1 = cdc_min_dedx + dedxbin1*cdc_dedx_step; + Float_t dedx2 = dedx1 + cdc_dedx_step; + + dedxcf = cf1 + (thisdedx - dedx1)*(cf2 - cf1)/(dedx2-dedx1); + + } else if (dedxbin1 == dedxbin2) { // interp theta only + + Float_t cf1 = CDC_DEDX_CORRECTION[dedxbin1][thetabin1]; + Float_t cf2 = CDC_DEDX_CORRECTION[dedxbin1][thetabin2]; + + Float_t theta1 = cdc_min_theta + thetabin1*cdc_theta_step; + Float_t theta2 = theta1 + cdc_theta_step; + + dedxcf = cf1 + (theta_deg - theta1)*(cf2 - cf1)/(theta2-theta1); + + } else { + + Float_t cf1 = CDC_DEDX_CORRECTION[dedxbin1][thetabin1]; + Float_t cf2 = CDC_DEDX_CORRECTION[dedxbin2][thetabin1]; + + Float_t dedx1 = cdc_min_dedx + dedxbin1*cdc_dedx_step; + Float_t dedx2 = dedx1 + cdc_dedx_step; + + Float_t cf3 = cf1 + (thisdedx - dedx1)*(cf2 - cf1)/(dedx2-dedx1); + + cf1 = CDC_DEDX_CORRECTION[dedxbin1][thetabin2]; + cf2 = CDC_DEDX_CORRECTION[dedxbin2][thetabin2]; + + dedx1 = cdc_min_dedx + dedxbin1*cdc_dedx_step; + dedx2 = dedx1 + cdc_dedx_step; + + Float_t cf4 = cf1 + (thisdedx - dedx1)*(cf2 - cf1)/(dedx2-dedx1); + + Float_t theta1 = cdc_min_theta + thetabin1*cdc_theta_step; + Float_t theta2 = theta1 + cdc_theta_step; + + dedxcf = cf3 + (theta_deg - theta1)*(cf4 - cf3)/(theta2-theta1); + + } + + return dedxcf; +} + + + +void checkmatrix(void){ + +Read_CDC_dEdx_CorrectionFile(); + +/* + // plot display options + gStyle->SetMarkerStyle(20); + gStyle->SetMarkerSize(0.5); + gStyle->SetOptFit(1); + gStyle->SetOptStat(0); + gStyle->SetOptLogz(1); +*/ + + // how many of each set of root files to read in. All of them would be 150 & 187. + int nchain1=150; //150; //number of files to add to the tchain (max 150) + int nchain2=187; //187; //187; + + // number of hits required per track + int minhits= 16; + + // first part of filenames + // char filestub1[] = "/raid12/gluex/rawdata2/Run050932/dedx_tree_nt_tree_050932"; + // char filestub2[] = "/raid12/gluex/rawdata2/Run050934/dedx_tree_nt_tree_050934"; + char filestub1[] = "/raid2/nsj/rawdata/Run050932/dedx_tree_4_18_1_nt_tree_050932"; + char filestub2[] = "/raid2/nsj/rawdata/Run050934/dedx_tree_4_18_1_nt_tree_50934"; + + // set up the chain of the root files which will be processed in turn + + TChain chain("dedx_tree"); + + for (int i=0; i78) continue; + if (r>1) continue; + + // find which histo to fill + + int ihisto = 0; + + while (theta > thetamax[ihisto]) ihisto++; + + // cout << "theta " << theta << " ihisto " << ihisto << " thetamax[ihisto] " << thetamax[ihisto] << endl; + + double new_dedx = dedx * (double)Correct_CDC_dEdx((float)theta, (float)dedx); + // fill the histograms + + if (ihisto < nhistos) { + if (charge > 0) { + dedxp_pos[ihisto]->Fill(p,new_dedx); + + dedxp_allpos->Fill(p,new_dedx); + dedxtheta_allpos->Fill(theta,new_dedx); + + } else { + dedxp_neg[ihisto]->Fill(p,new_dedx); + + dedxp_allneg->Fill(p,new_dedx); + dedxtheta_allneg->Fill(theta,new_dedx); + + } + } + + } + + + /* + + // Draw the histograms on the screen - VERY slow + + for (int i=30; i<31; i++) { //probably don't want 0 to nhistos unless nhistos is very small + + new TCanvas; + dedxp_pos[i]->Draw("colz"); + + new TCanvas; + dedxp_neg[i]->Draw("colz"); + + } + */ + + // create an output file & save the histograms into it + + TFile *f = new TFile(Form("matrix_applied_binned_histos_%ihits_%ifiles.root",minhits,nchain1+nchain2),"RECREATE"); + + for (int i=0; iWrite(); + for (int i=0; iWrite(); + dedxp_allpos->Write(); + dedxp_allneg->Write(); + dedxtheta_allpos->Write(); + dedxtheta_allneg->Write(); + // f->ls(); // list the file contents + f->Close(); + + +} diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/extendscalefactors.C b/src/plugins/Utilities/retired/dedx_tree/scripts/extendscalefactors.C new file mode 100644 index 0000000000..3ba5b9f01d --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/extendscalefactors.C @@ -0,0 +1,409 @@ + +void extendscalefactors(){ + + TFile *cffile = new TFile("dedx_correction_factors.root","READ"); + + TTree *t = (TTree*)gDirectory->Get("t"); + int nentries = (int)t->GetEntries(); + + int pid; + // 0: d 1: p 2: K 3: pi 4: combined + + const int THETA_MIN = 45; + const int THETA_MAX = 130; + int ntheta = 1+ THETA_MAX - THETA_MIN; + + const float DEDX_MIN=2.0; + const float DEDX_MAX=20.0; + + float p, theta; + float dedx, err_dedx, cf; + + t->SetBranchAddress("theta",&theta); + t->SetBranchAddress("dedx",&dedx); + t->SetBranchAddress("err_dedx",&err_dedx); + t->SetBranchAddress("f",&cf); + t->SetBranchAddress("p",&p); + t->SetBranchAddress("pid",&pid); + + + TTree *newt = new TTree("t","extended range of dedx correction factors"); + + newt->Branch("theta",&theta,"theta/F"); + newt->Branch("dedx",&dedx,"dedx/F"); + newt->Branch("err_dedx",&err_dedx,"err_dedx/F"); + newt->Branch("f",&cf,"f/F"); + // extras in case of problems + newt->Branch("p",&p,"p/F"); + newt->Branch("pid",&pid,"pid/I"); + + + // copy existing entries *****for pi and p only **** + + for (int i=0; iGetEntry(i); + if (pid==1||pid==3) newt->Fill(); + } + + + // the pions only run from 65 to 135 + // need to find the first and last pion cf for each p and copy it to every theta 45 - 64, 115-135 + // then fill in the dedx-min edge at 2. + + + // low dedx boundary + // scan up the column of pion dedx vs theta to find the cf for the lowest value of dedx + // ? if the lowest cf > 1.04 or has p > 2.6, leave a gap. + // ? The lowest cf band has a max at 1.04 w theta + + + + + + float lowestdedx[THETA_MAX+1]; + float cf_lowestdedx[THETA_MAX+1] = {0}; + + for (int i=THETA_MIN; i<= THETA_MAX; i++) lowestdedx[i] = 100.0; + + int itheta; + + for (int i=0; iGetEntry(i); + + if (pid != 3) continue; + if (dedx == 0) continue; + // if (cf>1.04) continue; + + itheta = int(theta); + // cout << dedx << " " << theta << " " << lowestdedx[itheta] << endl; + if (dedx < lowestdedx[itheta]) { + lowestdedx[itheta] = dedx; + cf_lowestdedx[itheta] = cf; + } + } + + + int firsttheta=0; + int lasttheta=0; + + + for (int itheta=THETA_MIN; itheta<=THETA_MAX; itheta++) { + + // cout << " min dedx at theta " << itheta << " degrees is " << lowestdedx[itheta] << endl; + + if (cf_lowestdedx[itheta] > 0) { + + if (firsttheta==0) firsttheta = itheta; + lasttheta = itheta; + + theta = float(itheta); + dedx = DEDX_MIN; + cf = cf_lowestdedx[itheta]; + pid = -1; + p = -1; + newt->Fill(); + cout << "filling in cf " << cf << " for dedx " << dedx << " theta " << theta << endl; + } + } + + + theta = (float)THETA_MIN; + cf = cf_lowestdedx[firsttheta]; + newt->Fill(); + + + theta = (float)THETA_MAX; + cf = cf_lowestdedx[lasttheta]; + newt->Fill(); + + + + // add high dedx boundary + + // make array of cf and dedx at each theta, then fit and extrapolate + + float allf[THETA_MAX+1][20]; + float alld[THETA_MAX+1][20]; + int alln[THETA_MAX+1] = {0}; + + float pionf[THETA_MAX+1][20]; + float piond[THETA_MAX+1][20]; + int pionn[THETA_MAX+1] = {0}; + + + for (int i=0; iGetEntry(i); + + int it = (int)theta; + + if (pid==1) { + alld[it][alln[it]] = dedx; + allf[it][alln[it]] = cf; + alln[it]++; + } else if (pid==3) { + piond[it][pionn[it]] = dedx; + pionf[it][pionn[it]] = cf; + pionn[it]++; + } + } + + TGraph *g; + TF1 *f = new TF1("f","pol1",0,16); + double fitstat; + + float lastgoodtheta, lastcfat20, lastcfat17, lastcfat15, lastcfat125,lastcfat4; // vals from highest theta before data run out + float lastgoodpitheta,firstcfat2,lastcfat2; + + + int counter = 0; + + for (int it = THETA_MIN; it<=THETA_MAX; it++) { + + if (alln[it]>1) { + + g = new TGraph(alln[it],alld[it],allf[it]); + + // if (counter==0) g->Draw("AP"); + + fitstat = g->Fit("f","Q"); + + if (!fitstat) { + + // evaluate the func for dedx=15, 20, and 12 if it is within the dip. + + pid = -1; + p=0; + theta = (float)it; + lastgoodtheta = theta; + + dedx = 20.0; + cf = f->Eval(dedx); + lastcfat20 = cf; + newt->Fill(); + + dedx = 17.5; + cf = f->Eval(dedx); + lastcfat17 = cf; + newt->Fill(); + + dedx = 15.0; + cf = f->Eval(dedx); + lastcfat15 = cf; + newt->Fill(); + + dedx = 12.5; //new + cf = f->Eval(dedx); + lastcfat125 = cf; + newt->Fill(); + + + dedx = 4.0; //new + cf = f->Eval(dedx); + lastcfat4 = cf; + newt->Fill(); + + + // if (theta<81) continue; + //if (theta>101) continue; + + if (theta >70 && theta < 112) { + dedx = 10.0; + cf = f->Eval(dedx); + newt->Fill(); + } + + + if (theta >80 && theta < 100) { + dedx = 9.0; + cf = f->Eval(dedx); + newt->Fill(); + } + + } // if !fitstat + + } + + /* + if (pionn[it]>1) { + + cout << "pion n " << pionn[it] << endl; + g = new TGraph(pionn[it],piond[it],pionf[it]); + + // if (counter==0) g->Draw("AP"); + + fitstat = g->Fit("f","0"); + + if (!fitstat) { + + // evaluate the func for dedx=2 + + pid = -1; + p=0; + theta = (float)it; + lastgoodpitheta = theta; + + dedx = 2.0; + cf = f->Eval(dedx); + if (firstcfat2==0) firstcfat2=cf; + lastcfat2 = cf; + newt->Fill(); + + } + + } // if alln[it]>0 elseif pin[it]>0 + */ + counter++; + if (counter==5) counter = 0; + + } + + + if (lastgoodtheta < THETA_MAX) { + + theta = THETA_MAX; + dedx = 20.0; + cf = lastcfat20; + newt->Fill(); + + dedx = 17.5; + cf = lastcfat15; + newt->Fill(); + + dedx =15.0; + cf = lastcfat15; + newt->Fill(); + + dedx =12.5; + cf = lastcfat125; + newt->Fill(); + + dedx =4.0; + cf = lastcfat4; + newt->Fill(); + + + } + + /* + theta = THETA_MIN; + dedx=2.0; + cf = firstcfat2; + newt->Fill(); + + theta = THETA_MAX; + cf = lastcfat2; + newt->Fill(); + */ + + //sides - first try without these. + + // how many mtm values are there? + + // 0.1 to 2.5 in steps of 0.05 + + float f1[5][280]={0}, t1[5][280]={90}, d1[5][280] = {0}; + float f2[5][280]={0}, t2[5][280]={90}, d2[5][280] = {0}; + + for (int j=0; j<5; j++) { + for (int i=0; i<280; i++) { + t1[j][i] = 90; + t2[j][i] = 90; + } + } + + int ip; //which array element to use for the p + + // find cf at first and last theta, for each mtm, protons & pions only + + for (int i=0; iGetEntry(i); + + ip = (int)(p*100); // so 0.35->35 + + if (theta < t1[pid][ip]) { + f1[pid][ip] = cf; + d1[pid][ip] = dedx; + t1[pid][ip] = theta; + } + + if (theta > t2[pid][ip]) { + f2[pid][ip] = cf; + d2[pid][ip] = dedx; + t2[pid][ip] = theta; + } + } + + // replicate them + + for (int ip=0; ip<280; ip++) { + + for (int thispid = 0; thispid<5; thispid++) { + + if (thispid == 2) continue; // no k+ + if (thispid == 0) continue; // no d + if (d1[thispid][ip] == 0) continue; + if (t1[thispid][ip] > 70) continue; // don't think there are any; too big to be useful + + // copy first value to theta_min degrees + pid = -1; + p = 0; + + + if (t1[thispid][ip] > THETA_MIN) { + + // if (thispid==3 && t1[thispid][ip]<=65) { + if (t1[thispid][ip]<=65) { + + theta = THETA_MIN; + cf = f1[thispid][ip]; + dedx = d1[thispid][ip]; + newt->Fill(); + + printf("for pid %i added cf %.1f from dedx %.1f p %.1f and theta %.1f at theta %.0f deg\n",thispid, cf, dedx, 0.01*ip, t1[thispid][ip],theta); + // cout << ip << " added " << dedx << " from " << t1[thispid][ip] << " at " << THETA_MIN << "deg "; + } + + } + + if (t2[thispid][ip] < THETA_MAX) { + + // if the last theta is < 100, duplicate cf from the first theta as they are mostly symmetrical + // cout << " last theta found is " << t2[thispid][ip] ; + + if ((t2[thispid][ip] >= 105 && thispid==3) || (t2[thispid][ip] >= 115 && thispid==1)) { + theta = THETA_MAX; + cf = f2[thispid][ip]; + dedx = d2[thispid][ip]; + + newt->Fill(); + + printf("\nfor pid %i added cf %.1f from dedx %.1f p %.1f and theta %.1f at theta %.0f deg\n",thispid, cf, dedx, 0.01*ip, t2[thispid][ip],theta); + // cout << " added " << dedx << " from " << t2[thispid][ip] << " degrees at " << THETA_MAX << "deg\n"; + + } else { + + /* + theta = THETA_MAX; + cf = f1[thispid][ip]; + dedx = d1[thispid][ip]; + + newt->Fill(); + + printf("for pid %i added cf %.1f from dedx %.1f p %.1f and theta %.1f at theta %.0f deg\n",thispid, cf, dedx, 0.01*ip, t1[thispid][ip],theta); + */ + + } + + } // end if (t2[thispid][ip] < THETA_MAX) { + + } // end for (int thispid = 0; thispid<5; thispid++) { + + } // end for (int ip=0; ip<280; ip++) { + + + TFile *newfile = new TFile("extended_dedx_correction_factors.root","RECREATE"); + newt->Write(); + newfile->Close(); + +} diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/fitbinnedhistoslg.C b/src/plugins/Utilities/retired/dedx_tree/scripts/fitbinnedhistoslg.C new file mode 100644 index 0000000000..7d95e1cc07 --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/fitbinnedhistoslg.C @@ -0,0 +1,564 @@ + +// Fits dedx vs p histograms in root file (name hard-coded) +// Uses a Gaussian to fit the proton peak and a Landau*Gaussian convolution for the pion peak +// Writes results into a tree and saves pictures of each fitted histo to simfitpics and simbadfitpics + +// uses these directories +// mkdir simfitpics +// mkdir simbadfitpics + +// mkdir simbadfitpics/p +// mkdir simbadfitpics/pi +// mkdir simfitpics/p +// mkdir simfitpics/pi + + +#include "langaus.C" // this contains the Landau*Gaussian fit function (obtained from root's example code ) + + +TF1 *onelangaufit(Char_t *FunName, TH1D *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF) +{ + + Int_t i; + + TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName); + if (ffitold) delete ffitold; + + TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4); + + ffit->SetNpx(100); + ffit->SetParameters(startvalues); + + ffit->SetParNames("Width","MP","Area","GSigma"); + + for (i=0; i<4; i++) { + ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]); + } + + his->Fit(FunName,"RB0WE"); // fit within specified range, use ParLimits, do not plot + + ffit->GetParameters(fitparams); // obtain fit parameters + for (i=0; i<4; i++) { + fiterrors[i] = ffit->GetParError(i); // obtain fit parameter errors + } + ChiSqr[0] = ffit->GetChisquare(); // obtain chi^2 + NDF[0] = ffit->GetNDF(); // obtain ndf + + return (ffit); // return fit function + +} + + +void fitbinnedhistoslg(void) { + + + // graph style options + gStyle->SetMarkerStyle(20); + gStyle->SetMarkerSize(0.5); + gStyle->SetOptFit(1111); + gStyle->SetOptStat(0); + gStyle->SetOptLogz(1); + + + + bool saveplots=1; // save plot png files + + if (saveplots) { + bool missingdir = 0; + char dirneeded[6][20] = {"simfitpics","simfitpics/p","simfitpics/pi","simbadfitpics","simbadfitpics/p","simbadfitpics/pi"}; + for (int i=0; i<6; i++) { + if (!gSystem->OpenDirectory(dirneeded[i])) { + printf("mkdir %s\n",dirneeded[i]); + missingdir = 1; + } + } + if (missingdir) gApplication->Terminate(); + } + + + bool sim=1; // if 1, fit simulated data - cuts are not theta-dependent + + // set to 1 to fit protons, pions, 0 to skip + bool fitp = 1; + bool fitpi = 1; + + + TFile *f; // file of binned histos + + if (sim) f = new TFile("/home/njarvis/online/dedx/16hits/sim_binned_histos_16hits_83files.root","READ"); + + if (!sim) f = new TFile("/home/njarvis/online/dedx/16hits/binned_histos_16hits_337files.root","READ"); + + + // 161 histos of each q+ and q- , dedxp+q+[n] + + // For 20+ hits, settled on range 45 degrees to 130 degrees. + // For 16+ hits, fitting 20 to 150 degrees + + const int nhistos=1+140-10; //1+130-45; //161; // # dedx vs p histos to process, ie # theta values, ultimately 161 + const int npbins=100; // # bin in p fitted, ie #p values , ultimately 100 + const int firsthisto=10; //35; //45 degrees //0; // ultimately 0, index of first histo to process + + + // p axis starts at -0.025, bin width is 0.05. + // bin 1 center is 0 + // bin 2 center is 0.05 + // bin 3 center is 0.1 + + float pstart = 0.05; // should be 0.05 // start bin value for momentum + float dp=0.05; // should be 0.05; // increment + + // TH2D *dedxp[nhistos]; + TH2D *hpos; + TH2D *hneg; + + + // fit functions + + TF1 *g = new TF1("gaus","gaus",0,25); + TF1 *lan = new TF1("landau","landau",0,25); + + g->SetLineColor(kBlue); + + TH1D *ypro; + + double fitstat1, fitstat2, pars1[3],pars2[3],chisq1,chisq2; + int ndof1,ndof2; + + float tstart,tstop; + + float theta; + float p; + + int pid; // 0: deuteron 1:proton 2:kaon 3:pion 4:merged p+pi band above 1.2GeV + + const float minpeakheight=10.0; // min acceptable peak height + + const int nneededtofit=150; // entries in a histo needed to fit. + const int nneededtoproject=500; // entries in a histo needed to cut. + + // error tolerances + + const float meanabserrtol=0.1; // keep fits with error on the mean less than this + const float meanfracerrtol=0.01; // keep fits with error on the mean less than this + + const float meanfracerrtol_d=0.1; // for d band: keep fits with error on the mean less than this + const float meanfracerrtol_k=0.1; // for k band: keep fits with error on the mean less than this + + // lines cutting between the bands. + // the parameter values depend on theta and are set later + + float pcutmax = 4.0; + TF1 *fp = new TF1("fp","exp([0]*(x + [1])) + [2]", 0.0, pcutmax); + + float kcutmax = 0.85; + TF1 *fk = new TF1("fk","exp([0]*(x + [1])) + [2]", 0.0, kcutmax); + + float picutmax = 0.425; + TF1 *fpi = new TF1("fpi","exp([0]*(x + [1])) + [2]", 0.0, picutmax); + + + // set up the output tree for the results + + TTree *t = new TTree("t","fit results"); + + t->Branch("p",&p,"p/F"); + t->Branch("theta",&theta,"theta/F"); + + double nentries; + double counts; // in the region for the pid of interest + t->Branch("n",&counts,"n/D"); + t->Branch("pid",&pid,"pid/I"); + + + // arrays to contain results from the fit function + double fitstat, fitpars[4],fiterrs[4]; + char parnames[4][10] ={"ht","mean","sigma","gsigma"}; + char errnames[4][10] ={"err_h","err_m","err_s","err_gs"}; + + + for (int i=0; i<4; i++) { + t->Branch(parnames[i],&fitpars[i],Form("%s/D",parnames[i])); + t->Branch(errnames[i],&fiterrs[i],Form("%s/D",errnames[i])); + } + + // LG fit params repeated temporarily + Double_t fpars_lg[4]; + Double_t fpe[4]; + + + for (int i=0; i<4; i++) t->Branch(Form("fp%i",i),&fpars_lg[i],Form("fp%i/D",i)); + for (int i=0; i<4; i++) t->Branch(Form("fpe%i",i),&fpe[i],Form("fpe%i/D",i)); + + + // text file to contain list of fit failures for the pi peak + // FILE *listbadpi = fopen("badpifits.txt","w"); + + // histogram bin number is theta - 10 + + for (int ihisto=firsthisto; ihistoSetParameters(-4.6, -0.98, 5.2); + fk->SetParameters(-4.6,-0.69,2.7); + fpi->SetParameters(-6.5, -0.43, 2.75); + + } else if (fabs(ihisto-80) < 10) { // ihisto 71-89 + fp->SetParameters(-4.7, -0.88, 4.0); // good for 90 degr + fk->SetParameters(-4.8,-0.65,2.3); // good for 90 degr + fpi->SetParameters(-6.5, -0.4, 2.4); + + } else if (fabs(ihisto-80) < 20) { // ihisto 61-70, 90-99 + fp->SetParameters(-4.6, -0.90, 4.8); // 80 degr + fk->SetParameters(-5.0,-0.64,3.0); + fpi->SetParameters(-6.5, -0.4, 2.8); + + } else if (fabs(ihisto-80) < 30) { // ihisto 51-60, 100-109 + fp->SetParameters(-4.6, -0.95, 4.8); // 70 degr + fk->SetParameters(-4.8,-0.69,2.7); + fpi->SetParameters(-6.5, -0.43, 2.75); + + } else if (fabs(ihisto-80) < 40) { // ihisto 41-50, 110-119 + fp->SetParameters(-4.6, -0.98, 5.2); // 51-60 degr + fk->SetParameters(-4.6,-0.7,2.3); // good at 60 degr + fpi->SetParameters(-6.5, -0.43, 2.75); + + } else if (fabs(ihisto-80) < 50) { // ihisto 31-40, 120-129 + fp->SetParameters(-4.6, -0.98, 5.2); // 41-50 degr + fk->SetParameters(-4.6,-0.69,2.7); // good at 50 degr + fpi->SetParameters(-6.5, -0.43, 2.75); + + } else { //if (fabs(ihisto-80) < 55) { // ihisto < 31, > 129 + fp->SetParameters(-4.6, -0.98, 5.2); // up to 40 degr + fk->SetParameters(-5.0, -0.7, 2.9); + fpi->SetParameters(-6.5, -0.43, 2.75); + + } + + // read in the q+ histogram from the root file of binned histos + + f->cd(); + + hpos = (TH2D*)gDirectory->Get(Form("dedxp_q+[%i]",ihisto)); + + cout << hpos->GetTitle() << " entries " << hpos->GetEntries() << endl; + + if (hpos->GetEntries() < nneededtoproject) continue; + + // new TCanvas; + // hpos->Draw("colz"); + + sscanf(hpos->GetTitle(),"theta %f to %f degrees",&tstart,&tstop); + theta = 0.5*(tstart+tstop); + + if (hpos->GetEntries() < nneededtoproject) continue; + + bool good; + bool checkfits = 0; // switch off fit quality checks + + // langaus fit params + Double_t fr[2]; + Double_t psv[4], pisv[4], pllo[4], plhi[4]; //, fpe[4]; //using variants of fp[] + // TF1* fitlg; + // Double_t fpars_lg[4] = {0}; + Double_t lgpeak, lgwidth; + Char_t funcname[100]; //name for current fit function + sprintf(funcname,"LanGaus"); + + Double_t chisqr; + int ndf; + + // scan through p bins + for (int ip = 0; ip < npbins; ip++) { + + p = pstart + ip*dp; + + double pbin = hpos->GetXaxis()->FindBin(p); + ypro = hpos->ProjectionY("p",pbin,pbin); + + ypro->SetTitle(Form("q+ %s momentum %.2f GeV/c",hpos->GetTitle(),p)); + + nentries = ypro->GetEntries(); + + //cout << "p " << p << " theta " << theta << " " << "entries " << nentries << endl; + + if (nentries < nneededtofit) continue; + + double yrangemax = ypro->GetYaxis()->GetXmax(); + + // use dE/dx cut func to separate the bands, find the max bin either side + + new TCanvas; ypro->DrawCopy(); + + // cut functions + + double pitop = fpi->Eval(p,0,0); + double ktop = fk->Eval(p,0,0); + double ptop = fp->Eval(p,0,0); + double dtop = ptop + 6; + + if (p > picutmax) pitop = ktop; + + + + // find counts between these bins + + double pcounts = ypro->Integral(ypro->GetXaxis()->FindBin(ktop),ypro->GetXaxis()->FindBin(ptop)); + double picounts = ypro->Integral(0,ypro->GetXaxis()->FindBin(pitop)); + + + // set fit results to 0 in case the fit is switched off or there aren't enough counts + fitstat = 1; + + for (int i=0; i<4; i++) { + fitpars[i] = 0; + fiterrs[i] = 0; + fpars_lg[i] = 0; + fpe[i] = 0; + } + + + + // ------------------ protons ------------------------------- + + printf("\n\n\n p %.2f GeV/c theta %2f degrees\n\n",p,theta); + printf("Counts in peak regions: p %f pi %f\n", pcounts, picounts); + + bool unfittable = 0; + + // unfittable because the p and pi are too close + if (fabs(p-0.75)<0.01 && (theta > 67 && theta < 82)) unfittable = 1; + if (fabs(p-0.80)<0.01 && (theta > 66 && theta < 83)) unfittable = 1; + if (fabs(p-0.85)<0.01 && (theta > 64 && theta < 80)) unfittable = 1; + + // unfittable due to stats + if (fabs(p-0.80)<0.01 && (theta > 99 && theta < 102)) unfittable = 1; + if (fabs(p-0.85)<0.01 && (theta > 79 && theta < 84)) unfittable = 1; + + + pid = 1; + + if (fitp && !unfittable && p<=kcutmax && pcounts>nneededtofit) { + + // find ktop and ptop, fit data between ktop and ptop with gauss + + g->SetParameter(1, 0.5*(ktop+ptop)); + g->SetParameter(2,1); + + fitstat = ypro->Fit(g,"","",ktop,ptop); + + //cout << "p fit, ktop " << ktop << " ptop " << ptop << endl; + + good = 1; + + if (!fitstat) { + + g->GetParameters(&fitpars[0]); + for (int i=0; i<3; i++) fiterrs[i] = g->GetParError(i); + + // check peak is within range + if (fitpars[1] < ktop) good = 0; + if (fitpars[1] > ptop) good = 0; + + // check peak height + if (fitpars[0] < minpeakheight ) good = 0; //not enough counts + + // check error:param ratio + for (int i=0; i<3; i++) if (fiterrs[i] >= fitpars[i]) good = 0; + + + if (good && checkfits) { + + if (fitpars[1] < ktop + 1.0*fitpars[2] ) good = 0; // peak & cut within 1sigma + if (fitpars[1] < ktop + 1.0*fitpars[2] ) printf("mean - ktop %.2f < limit %.2f\n",ktop,1.0*fitpars[2]); + + if (fiterrs[1] >= meanabserrtol) printf("error on the mean %.2f > limit %.2f\n",fiterrs[1],meanabserrtol); + if (fiterrs[1] >= meanabserrtol) good = 0; + + if (fiterrs[1]/fitpars[1] >= meanfracerrtol) printf("err:mean %.2f > limit %.2f\n",fiterrs[1]/fitpars[1],meanfracerrtol); + if (fiterrs[1]/fitpars[1] >= meanfracerrtol) good = 0; + } + } else { + + good = 0; + + } + + if (saveplots) ypro->GetYaxis()->SetRangeUser(0,1.5*fitpars[0]); + if (saveplots) gPad->Update(); + + if (good) { + + counts = pcounts; + t->Fill(); + + if (saveplots) gPad->SaveAs(Form("simfitpics/p/dedx_p%.2f_theta%.0f.png",p,theta)); + //if (saveplots) gPad->SaveAs(Form("simfitpics/p/theta%.0f_dedx_p%.2f.png",theta,p)); + } else { + if (saveplots) gPad->SaveAs(Form("simbadfitpics/p/dedx_p%.2f_theta%.0f.png",p,theta)); + //if (saveplots) gPad->SaveAs(Form("simbadfitpics/p/theta%.0f_dedx_p%.2f.png",theta,p)); + } + + } + + + + // ------------------ pions ------------------------------- + + pid = 3; + + unfittable = 0; + + // unfittable because the p and pi are too close + if (fabs(p-0.75)<0.01 && (theta > 85 && theta < 90)) unfittable = 1; + //if (fabs(p-0.70)<0.01 && (theta > 87)) unfittable = 1; + if (fabs(p-0.80)<0.01 && (theta > 78 && theta < 85)) unfittable = 1; + if (fabs(p-0.85)<0.01 && (theta<63 || theta>75)) unfittable = 1; + + // unfittable due to stats + // if (fabs(p-0.55)<0.01 && (theta > 117 && theta < 119)) unfittable = 1; + // if (fabs(p-0.6)<0.01 && (theta > 105 && theta < 108)) unfittable = 1; + + + if (fitpi && !unfittable && (p <= kcutmax ) && picounts>nneededtofit) { + + // printf("\n\n\n p %.2f GeV/c theta %2f degrees\n\n",p,theta); + // printf("Counts in peak regions: d: %f p %f K %f pi %f\n", dcounts, pcounts, kcounts, picounts); + // find ktop and ptop, fit data between ktop and ptop with gauss + + g->SetParameter(1, 0.5*pitop); + g->SetParameter(2,1); + fitstat = ypro->Fit(g,"","",1,pitop); + g->GetParameters(&fitpars[0]); + + if (!fitstat && (fitpars[1]GetRMS(); //lpars[2]; //1.8; + pisv[1] = 0.9*fitpars[1]; + pisv[2] = 2.0*fitpars[0]; //0.1*picounts; + pisv[3] = fitpars[2]; //0.4; + } else { + pisv[0] = 0.2; //0.1*fitpars[2]; //0.1*dneg->GetRMS(); //lpars[2]; //1.8; + pisv[1] = 0.5*(1.0+pitop); + pisv[2] = 0.1*picounts; + pisv[3] = 0.4; + } + + + fr[0] = 1; + fr[1]=pitop; //3.0*dneg->GetMean(); + + good = 1; + + // TF1 *langaufit(TH1F *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF) + // par[0]=Width (scale) parameter of Landau density + // par[1]=Most Probable (MP, location) parameter of Landau density + // par[2]=Total area (integral -inf to inf, normalization constant) + // par[3]=Width (sigma) of convoluted Gaussian function + + cout << "initial guesses: " << pisv[0] << " " << pisv[1] << " " << pisv[2] << " " << pisv[3] << endl; + + pllo[0]=0.05; pllo[1]=1.0; pllo[2]=pisv[2]/5.; pllo[3]=0.05; //pisv[3]*0.66; + plhi[0]=1.0; plhi[1]=pitop; plhi[2]=pisv[2]*5.; plhi[3]=1.0; //pisv[3]*1.33; + + TF1* fitlg = onelangaufit(funcname,ypro,fr,pisv,pllo,plhi,fpars_lg,fpe,&chisqr,&ndf); + + fitlg->Draw("same"); + + printf ("Chi sq %f NDF %i Chi sq/NDF %f\n",chisqr,ndf,chisqr/(float)ndf); + + // check peak is inside fit range + if (fitpars[1] < fr[0]) good = 0; + if (fitpars[1] > fr[1]) good = 0; + + // check peak height + fitpars[0] = fitlg->Eval(fitpars[1]); + if (fitpars[0]= fpars_lg[j]) good = 0; + + + + if (good && checkfits) { + + for (int j=0; j<4; j++) if (fpe[j] > 0.2*fpars_lg[j]) good = 0; + for (int j=0; j<4; j++) if (fpe[j] > 0.2*fpars_lg[j]) printf("*** param %i error > 0.1 x param val ***\n",j); + + + if (fpe[0] >= 2.0*meanabserrtol) good = 0; + if (fpe[0] >= 2.0*meanabserrtol) printf("*** Width error %.2f is too large ***\n",fpe[0]); + + if (fpe[1] >= meanabserrtol) good = 0; + if (fpe[1] >= meanabserrtol) printf("*** Peak error %.2f is too large ***\n",fpe[1]); + } + + if (good) { + + counts = picounts; + + // fitpars[0] = fpars_lg[2]; // par[2]=Total area (integral -inf to inf, normalization constant) + fitpars[1] = fpars_lg[1]; //lgpeak; //fpars_lg[0]; + fitpars[2] = 0.5*lgwidth; //fpars_lg[0]; + fitpars[3] = 0.5*fpars_lg[3]; // half gaussian sigma + + fiterrs[0] = 0; //fpe[2]; // total area + fiterrs[1] = fpe[1]; // mpv error // par[1]= MPV of Landau density + fiterrs[2] = 0.5*fpe[0]; // landau width // par[0]=Width (scale) param of Landau density + fiterrs[3] = 0.5*fpe[3]; // par[3]=Width (sigma) of convoluted Gaussian function + + } + + + if (saveplots) ypro->GetYaxis()->SetRangeUser(0,1.0*fpars_lg[2]); + if (saveplots) ypro->GetXaxis()->SetRangeUser(0,10); + if (saveplots) gPad->Update(); + + if (good) { + + t->Fill(); + + if (saveplots) gPad->SaveAs(Form("simfitpics/pi/dedx_p%.2f_theta%.0f.png",p,theta)); + //if (saveplots) gPad->SaveAs(Form("simfitpics/pi/theta%.0f_dedx_p%.2f.png",theta,p)); + + } else { + + // fprintf(listbadpi,"theta %.0f p %.2f \n",theta,p); + + if (saveplots) gPad->SaveAs(Form("simbadfitpics/pi/dedx_p%.2f_theta%.0f.png",p,theta)); + //if (saveplots) gPad->SaveAs(Form("simbadfitpics/pi/theta%.0f_dedx_p%.2f.png",theta,p)); + } + + } + + + } // ip loop (momentum) + + } // ihisto loop (theta) + + + // save the fit results to a new tree file + + TFile *outfile; + + if (sim) outfile = new TFile("fittedsimresults.root","RECREATE"); + if (!sim) outfile = new TFile("fitresults.root","RECREATE"); + outfile->cd(); + t->Write(); + + outfile->Close(); + + // fclose(listbadpi); + +} + + diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/langaus.C b/src/plugins/Utilities/retired/dedx_tree/scripts/langaus.C new file mode 100644 index 0000000000..990b0dd573 --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/langaus.C @@ -0,0 +1,273 @@ +//----------------------------------------------------------------------- +// +// Convoluted Landau and Gaussian Fitting Function +// (using ROOT's Landau and Gauss functions) +// +// Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at) +// Adapted for C++/ROOT by H.Pernegger (Heinz.Pernegger@cern.ch) and +// Markus Friedl (Markus.Friedl@cern.ch) +// +// to execute this example, do: +// root > .x langaus.C +// or +// root > .x langaus.C++ +// +//----------------------------------------------------------------------- + +#include "TH1.h" +#include "TF1.h" +#include "TROOT.h" +#include "TStyle.h" +#include "TMath.h" + +Double_t langaufun(Double_t *x, Double_t *par) { + + //Fit parameters: + //par[0]=Width (scale) parameter of Landau density + //par[1]=Most Probable (MP, location) parameter of Landau density + //par[2]=Total area (integral -inf to inf, normalization constant) + //par[3]=Width (sigma) of convoluted Gaussian function + // + //In the Landau distribution (represented by the CERNLIB approximation), + //the maximum is located at x=-0.22278298 with the location parameter=0. + //This shift is corrected within this function, so that the actual + //maximum is identical to the MP parameter. + + // Numeric constants + Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2) + Double_t mpshift = -0.22278298; // Landau maximum location + + // Control constants + Double_t np = 100.0; // number of convolution steps + Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas + + // Variables + Double_t xx; + Double_t mpc; + Double_t fland; + Double_t sum = 0.0; + Double_t xlow,xupp; + Double_t step; + Double_t i; + + + // MP shift correction + mpc = par[1] - mpshift * par[0]; + + // Range of convolution integral + xlow = x[0] - sc * par[3]; + xupp = x[0] + sc * par[3]; + + step = (xupp-xlow) / np; + + // Convolution integral of Landau and Gaussian by sum + for(i=1.0; i<=np/2; i++) { + xx = xlow + (i-.5) * step; + fland = TMath::Landau(xx,mpc,par[0]) / par[0]; + sum += fland * TMath::Gaus(x[0],xx,par[3]); + + xx = xupp - (i-.5) * step; + fland = TMath::Landau(xx,mpc,par[0]) / par[0]; + sum += fland * TMath::Gaus(x[0],xx,par[3]); + } + + return (par[2] * step * sum * invsq2pi / par[3]); +} + + + +TF1 *langaufit(TH1F *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF) +{ + // Once again, here are the Landau * Gaussian parameters: + // par[0]=Width (scale) parameter of Landau density + // par[1]=Most Probable (MP, location) parameter of Landau density + // par[2]=Total area (integral -inf to inf, normalization constant) + // par[3]=Width (sigma) of convoluted Gaussian function + // + // Variables for langaufit call: + // his histogram to fit + // fitrange[2] lo and hi boundaries of fit range + // startvalues[4] reasonable start values for the fit + // parlimitslo[4] lower parameter limits + // parlimitshi[4] upper parameter limits + // fitparams[4] returns the final fit parameters + // fiterrors[4] returns the final fit errors + // ChiSqr returns the chi square + // NDF returns ndf + + Int_t i; + Char_t FunName[100]; + + sprintf(FunName,"Fitfcn_%s",his->GetName()); + + TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName); + if (ffitold) delete ffitold; + + TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4); + ffit->SetParameters(startvalues); + ffit->SetParNames("Width","MP","Area","GSigma"); + + for (i=0; i<4; i++) { + ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]); + } + + his->Fit(FunName,"RB0"); // fit within specified range, use ParLimits, do not plot + + ffit->GetParameters(fitparams); // obtain fit parameters + for (i=0; i<4; i++) { + fiterrors[i] = ffit->GetParError(i); // obtain fit parameter errors + } + ChiSqr[0] = ffit->GetChisquare(); // obtain chi^2 + NDF[0] = ffit->GetNDF(); // obtain ndf + + return (ffit); // return fit function + +} + + +Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) { + + // Seaches for the location (x value) at the maximum of the + // Landau-Gaussian convolute and its full width at half-maximum. + // + // The search is probably not very efficient, but it's a first try. + + Double_t p,x,fy,fxr,fxl; + Double_t step; + Double_t l,lold; + Int_t i = 0; + Int_t MAXCALLS = 10000; + + + // Search for maximum + + p = params[1] - 0.1 * params[0]; + step = 0.05 * params[0]; + lold = -2.0; + l = -1.0; + + + while ( (l != lold) && (i < MAXCALLS) ) { + i++; + + lold = l; + x = p + step; + l = langaufun(&x,params); + + if (l < lold) + step = -step/10; + + p += step; + } + + if (i == MAXCALLS) + return (-1); + + maxx = x; + + fy = l/2; + + + // Search for right x location of fy + + p = maxx + params[0]; + step = params[0]; + lold = -2.0; + l = -1e300; + i = 0; + + + while ( (l != lold) && (i < MAXCALLS) ) { + i++; + + lold = l; + x = p + step; + l = TMath::Abs(langaufun(&x,params) - fy); + + if (l > lold) + step = -step/10; + + p += step; + } + + if (i == MAXCALLS) + return (-2); + + fxr = x; + + + // Search for left x location of fy + + p = maxx - 0.5 * params[0]; + step = -params[0]; + lold = -2.0; + l = -1e300; + i = 0; + + while ( (l != lold) && (i < MAXCALLS) ) { + i++; + + lold = l; + x = p + step; + l = TMath::Abs(langaufun(&x,params) - fy); + + if (l > lold) + step = -step/10; + + p += step; + } + + if (i == MAXCALLS) + return (-3); + + + fxl = x; + + FWHM = fxr - fxl; + return (0); +} + +void langaus() { + // Fill Histogram + Int_t data[100] = {0,0,0,0,0,0,2,6,11,18,18,55,90,141,255,323,454,563,681, + 737,821,796,832,720,637,558,519,460,357,291,279,241,212, + 153,164,139,106,95,91,76,80,80,59,58,51,30,49,23,35,28,23, + 22,27,27,24,20,16,17,14,20,12,12,13,10,17,7,6,12,6,12,4, + 9,9,10,3,4,5,2,4,1,5,5,1,7,1,6,3,3,3,4,5,4,4,2,2,7,2,4}; + TH1F *hSNR = new TH1F("snr","Signal-to-noise",400,0,400); + + for (Int_t i=0; i<100; i++) hSNR->Fill(i,data[i]); + + // Fitting SNR histo + printf("Fitting...\n"); + + // Setting fit range and start values + Double_t fr[2]; + Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4]; + fr[0]=0.3*hSNR->GetMean(); + fr[1]=3.0*hSNR->GetMean(); + + pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.4; + plhi[0]=5.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=5.0; + sv[0]=1.8; sv[1]=20.0; sv[2]=50000.0; sv[3]=3.0; + + Double_t chisqr; + Int_t ndf; + TF1 *fitsnr = langaufit(hSNR,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf); + + Double_t SNRPeak, SNRFWHM; + langaupro(fp,SNRPeak,SNRFWHM); + + printf("Fitting done\nPlotting results...\n"); + + // Global style settings + gStyle->SetOptStat(1111); + gStyle->SetOptFit(111); + gStyle->SetLabelSize(0.03,"x"); + gStyle->SetLabelSize(0.03,"y"); + + hSNR->GetXaxis()->SetRange(0,70); + hSNR->Draw(); + fitsnr->Draw("lsame"); +} + diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/make_histo_into_text_file.C b/src/plugins/Utilities/retired/dedx_tree/scripts/make_histo_into_text_file.C new file mode 100644 index 0000000000..fbef1dd450 --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/make_histo_into_text_file.C @@ -0,0 +1,210 @@ +{ + + TTree* t = (TTree*)gDirectory->Get("t"); + + int nentries = (int)t->GetEntries(); + float theta, dedx, f, p; + int pid; + t->SetBranchAddress("dedx",&dedx); + t->SetBranchAddress("theta",&theta); + t->SetBranchAddress("f",&f); + t->SetBranchAddress("pid",&pid); + t->SetBranchAddress("p",&p); + + const int THETA_MIN = 21; + const int THETA_MAX = 136; + int ntheta = 1 + THETA_MAX - THETA_MIN ; + + // theta 22 to 135, extended to 21 to 136 + // dedx 2.30382 to 15.0024, extended to 20 + // histo dedx 2.30740 to 20 + + + TCanvas *c = new TCanvas("c","Graph2D example",0,0,800,600); + Double_t x, y, z; + // Int_t np = 200; + + const int npy = 1+ 100; + + TGraph2D *gtot = new TGraph2D(); + gtot->SetNpx(ntheta); + gtot->SetNpy(npy); + + TGraph2D *gprotons = new TGraph2D(); + gprotons->SetNpx(ntheta); + gprotons->SetNpy(npy); + + TGraph2D *gpions = new TGraph2D(); + gpions->SetNpx(ntheta); + gpions->SetNpy(npy); + + TGraph2D *gdeuterons = new TGraph2D(); + gdeuterons->SetNpx(ntheta); + gdeuterons->SetNpy(npy); + + TGraph2D *gcombo = new TGraph2D(); + gcombo->SetNpx(ntheta); + gcombo->SetNpy(npy); + + TGraph2D *gfakes = new TGraph2D(); + gfakes->SetNpx(ntheta); + gfakes->SetNpy(npy); + + int np=0; + int nf=0; + int nd=0; + int npi=0; + int nc=0; + + int nt=0; // use for drawing more than one particle, just draw over one in a diff colour. + // to get the axis range right + + float dcut = 0.46; + + for (Int_t N=0; NGetEntry(N); + x = theta; //2*P*(r->Rndm(N))-P; + y = dedx; //2*P*(r->Rndm(N))-P; + z = f; //(sin(x)/x)*(sin(y)/y)+0.2; + if (pid==-1) gfakes->SetPoint(nf,x,y,z); + if (pid==0 && p < dcut) gdeuterons->SetPoint(nd,x,y,z); + if (pid==1) gprotons->SetPoint(np,x,y,z); + if (pid==3) gpions->SetPoint(npi,x,y,z); + if (pid==4) gcombo->SetPoint(nc,x,y,z); + + if (fabs(pid)!=0 && pid!=4) gtot->SetPoint(nt,x,y,z); + if (fabs(pid) != 0 && pid !=4) nt++; + + // if (pid==0 && pSetPoint(nt,x,y,z); + // if (pid==0 && pSetPalette(1); + gtot->SetTitle("dE/dx correction factor;theta;dedx;f"); + gtot->SetMarkerColor(0); + gtot->Draw("surf1"); + + TH2D* histo = gtot->GetHistogram(); + + new TCanvas; + histo->Draw("surf1"); + + cout << npy << endl; + + char filename[100]; + sprintf(filename,"surface_histo_%i_%i.png",ntheta,npy); + + c1->SaveAs(filename); + + sprintf(filename,"surface_histo_%i_%i.root",ntheta,npy); + + TFile* newfile = new TFile(filename,"RECREATE"); + histo->Write(); + newfile->Close(); + + // histo has npy bins in y and ntheta bins in x + + // bins are ordered underflow data overflow + // ie 2 extra bins per histo + + double thetamin = histo->GetXaxis()->GetXmin(); + double thetamax = histo->GetXaxis()->GetXmax(); + + double dedxmin = histo->GetYaxis()->GetXmin(); + double dedxmax = histo->GetYaxis()->GetXmax(); + + double thetastep = (thetamax - thetamin) / (float)(ntheta - 1); + double dedxstep = (dedxmax - dedxmin) / (float)(npy - 1); + + sprintf(filename,"dedx_corrections_%i_%i.txt",ntheta,npy); + FILE *outfile = fopen(filename,"w"); + fprintf(outfile,"%i values of theta\n",ntheta); //histo->GetNbinsX()); + fprintf(outfile,"%f min theta\n",thetamin); + fprintf(outfile,"%f max theta\n",thetamax); + fprintf(outfile,"%f theta step\n",thetastep); + + fprintf(outfile,"%i values of dedx\n",npy); + fprintf(outfile,"%f min dedx\n",dedxmin); + fprintf(outfile,"%f max dedx\n",dedxmax); + fprintf(outfile,"%f dedx step\n",dedxstep); + + fprintf(outfile,"\n"); + + + int totalbins = (ntheta+2)*(npy+2); + int ix = 0; + + for (int ibin = ntheta+2; ibin0 && ix<=ntheta) fprintf(outfile,"%f\n",histo->GetBinContent(ibin)); + ix++; + if (ix== ntheta+2) ix=0; + } + + fclose(outfile); + + /* + int ibin=0; + + cout << histo->GetNbinsX() << "x bins" << endl; + cout << histo->GetNbinsY() << "y bins" << endl; + + cout << endl; + + cout << "first bin 0\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + ibin=ntheta-1; + cout << "last bin in bottom row " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + + + ibin = ntheta; + cout << "bin ntheta " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + ibin = ntheta+1; + cout << "bin ntheta+1: " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + ibin = ntheta+2; + cout << "bin ntheta+2: " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + + ibin = ntheta+3; + cout << "bin ntheta+3: " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + + ibin = ntheta+4; + cout << "bin ntheta+4: " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + + // cout << histo->GetXaxis()->GetBinLowEdge(ibin); + // cout << histo->GetYaxis()->GetBinLowEdge(ibin); + // cout << histo->GetBinContent(ibin); + + + ibin = (ntheta+0)*npy; + cout << "last bin " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + + ibin = (ntheta+1)*npy; + cout << "last overflow bin " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + + +*/ + +} diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/makebinnedhistos.C b/src/plugins/Utilities/retired/dedx_tree/scripts/makebinnedhistos.C new file mode 100644 index 0000000000..56efbc718a --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/makebinnedhistos.C @@ -0,0 +1,158 @@ + +{ + + // plot display options + gStyle->SetMarkerStyle(20); + gStyle->SetMarkerSize(0.5); + gStyle->SetOptFit(1); + gStyle->SetOptStat(0); + gStyle->SetOptLogz(1); + + // how many of each set of root files to read in. All of them would be 150 & 187. + int nchain1=1; //150; //number of files to add to the tchain (max 150) + int nchain2=0; //187; + + // number of hits required per track + int minhits= 16; + + // first part of filenames + // char filestub1[] = "/raid12/gluex/rawdata2/Run050932/dedx_tree_nt_tree_050932"; + // char filestub2[] = "/raid12/gluex/rawdata2/Run050934/dedx_tree_nt_tree_050934"; + char filestub1[] = "/raid2/nsj/rawdata/Run050932/dedx_tree_4_18_1_nt_tree_050932"; + char filestub2[] = "/raid2/nsj/rawdata/Run050934/dedx_tree_4_18_1_nt_tree_50934"; + + // set up the chain of the root files which will be processed in turn + + TChain chain("dedx_tree"); + + for (int i=0; i78) continue; + if (r>1) continue; + + // find which histo to fill + + int ihisto = 0; + + while (theta > thetamax[ihisto]) ihisto++; + + // cout << "theta " << theta << " ihisto " << ihisto << " thetamax[ihisto] " << thetamax[ihisto] << endl; + + + // fill the histograms + + if (ihisto < nhistos) { + if (charge > 0) { + dedxp_pos[ihisto]->Fill(p,dedx); + } else { + dedxp_neg[ihisto]->Fill(p,dedx); + } + } + + } + + + /* + + // Draw the histograms on the screen - VERY slow + + for (int i=30; i<31; i++) { //probably don't want 0 to nhistos unless nhistos is very small + + new TCanvas; + dedxp_pos[i]->Draw("colz"); + + new TCanvas; + dedxp_neg[i]->Draw("colz"); + + } + */ + + // create an output file & save the histograms into it + + TFile *f = new TFile(Form("binned_histos_%ihits_%ifiles.root",minhits,nchain1+nchain2),"RECREATE"); + + for (int i=0; iWrite(); + for (int i=0; iWrite(); + // f->ls(); // list the file contents + f->Close(); + + +} diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/makescalefactors.C b/src/plugins/Utilities/retired/dedx_tree/scripts/makescalefactors.C new file mode 100644 index 0000000000..91a63b7b45 --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/makescalefactors.C @@ -0,0 +1,218 @@ + +void makescalefactors(){ + + // sim file has 12 entries for p, 16 for pi+ + // read these into an array + // then read in the real file, make new tree + // dedx, theta, correction-factor,pid + // include pid for tracking down weird points + + // also check to make sure that there aren't enough d in sim to fit + + + TFile *f = new TFile("fittedsimresults.root","READ"); + TTree *simt = (TTree*)gDirectory->Get("t"); + int simentries = (int)simt->GetEntries(); + + const int simarraysize=1000; //992 proton entries. only using 90 degrees for sim though. + + int pid; + // 0: d 1: p 2: K 3: pi 4: combined + + float p, theta; + double mean, err, ht; + + simt->SetBranchAddress("pid",&pid); + simt->SetBranchAddress("p",&p); + simt->SetBranchAddress("ht",&ht); + simt->SetBranchAddress("theta",&theta); + simt->SetBranchAddress("mean",&mean); + simt->SetBranchAddress("err_m",&err); + + float sim_proton_theta[simarraysize]; + float sim_pion_theta[simarraysize]; + + float sim_proton_mtm[simarraysize]; + float sim_pion_mtm[simarraysize]; + + double sim_proton_dedx[simarraysize]; + double sim_pion_dedx[simarraysize]; + + int np=0; + int npi=0; + + bool usethis; + + // load in sim data for limited range of theta + + for (int i=0; iGetEntry(i); + + // round p to 2 sig figs to make code easier + p = 0.01*int(p*100); + + if (pid==1 && (int)theta>44 && (int)theta<131 && ht>25) { + sim_proton_mtm[np] = p; + sim_proton_dedx[np] = mean; + sim_proton_theta[np] = theta; + np++; + + } else if (pid==3 && (int)theta==90 && ht>100) { + + sim_pion_mtm[npi] = p; + sim_pion_dedx[npi] = mean; + sim_pion_theta[npi] = theta; + npi++; + + } + } + + f->Close(); + + + + f = new TFile("fitresults.root","READ"); + TTree *t = (TTree*)gDirectory->Get("t"); + int nentries = (int)t->GetEntries(); + + t->SetBranchAddress("ht",&ht); + t->SetBranchAddress("pid",&pid); + t->SetBranchAddress("p",&p); + t->SetBranchAddress("theta",&theta); + t->SetBranchAddress("mean",&mean); + t->SetBranchAddress("err_m",&err); + + + + TTree *newt = new TTree("t","dedx correction factors"); + + float dedx, err_dedx, cf; + + newt->Branch("theta",&theta,"theta/F"); + newt->Branch("dedx",&dedx,"dedx/F"); + newt->Branch("err_dedx",&err_dedx,"err_dedx/F"); + newt->Branch("f",&cf,"f/F"); + // extras in case of problems + newt->Branch("p",&p,"p/F"); + newt->Branch("pid",&pid,"pid/I"); + + double firstdedx; // use this to save proton dedx at 45 or 67 degrees + double lastdedx; + double lookfortheta; + bool found; + + int itheta; + float simp; + + for (int i=0; iGetEntry(i); + + // round p to 2 sig figs to make code easier + p = 0.01*int(p*100); + itheta = (int) theta; + + found = 0; + usethis = 1; + + // printf("looking for pid %i and mtm %f\n",pid,p); + + // if (pid == 1 && itheta == 66 && p ==0.8) continue; // remove outlier + + if (pid ==1 && itheta >= 45 && itheta <= 130 && p > 0.349) { + + if (ht < 60) usethis = 0; + + if (p==0.35 && (itheta < 67)) usethis = 0; + if (p==0.35 && (itheta > 83)) usethis = 0; + + if (usethis) { + + found = 0; + firstdedx = 0; + + for (int j = 0; j < np; j++) { // find dedx for the first theta for the current p + + if (p != sim_proton_mtm[j]) continue; + + lastdedx = sim_proton_dedx[j]; // they are ordered by theta and then by p + + if (p == 0.35 && (int)sim_proton_theta[j]==67) { + firstdedx = sim_proton_dedx[j]; + } else if ((int)sim_proton_theta[j]==45) { + firstdedx = sim_proton_dedx[j]; + } + + } + + printf("set firstdedx to %f for p %f \n",firstdedx,p); + + for (int j = 0; j < np; j++) { + + if (p != sim_proton_mtm[j]) continue; + + if (itheta == (int)sim_proton_theta[j]) { + found = 1; + dedx = mean; + err_dedx = err; + cf = sim_proton_dedx[j]/dedx; + } + if (found) printf("Found a direct match for p %.2f and theta %.0f\n",p,theta); + if (found) break; + } + + if (!found) { + found = 1; + dedx = mean; + err_dedx = err; + cf = lastdedx/dedx; //always present + printf("Using lastdedx %f for p %.2f and theta %.0f, firstdedx is %f\n",lastdedx,p,theta,firstdedx); + } + + if (!found) printf("Could not find a sim dedx for p %.2f and theta %.0f\n",p,theta); + + } + + } else if (pid == 3 && itheta >= 60 && itheta <=110 && p >= 0.15 && p <= 0.45) { + + if (p==0.15 && (itheta < 70)) usethis = 0; + + if (ht < 100) usethis = 0; + + if (usethis) { + + found = 0; + + lookfortheta = 90; + + for (int j = 0; j < npi; j++) { + if ((p == sim_pion_mtm[j]) && (sim_pion_theta[j] ==lookfortheta)) { + found = 1; + dedx = mean; + err_dedx = err; + cf = sim_pion_dedx[j]/dedx; + } + if (found) break; + } + + if (!found) printf("Did not find mtm to match %.2f and theta %.0f for pid %i\n",p,theta,pid); + } + + } + + if (found) newt->Fill(); + + // read in from t + // search arrays for matching p and pid + // calc correction factor f + // fill new tree + // write out to new file + } + + + TFile *newfile = new TFile("dedx_correction_factors.root","RECREATE"); + newt->Write(); + newfile->Close(); + +} diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/makesimsinglebinnedhistos.C b/src/plugins/Utilities/retired/dedx_tree/scripts/makesimsinglebinnedhistos.C new file mode 100644 index 0000000000..b7266bcfca --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/makesimsinglebinnedhistos.C @@ -0,0 +1,175 @@ + +{ + + + // how many of each set of root files to read in. + int nchain=5; // number of files to add to the tchain + int nchain2=9; + int nchain3=8; + int nchain4=8; + int nchain5=8+ 53 - 8; //up to _083 + + int totalnchain = nchain + nchain2+nchain3+nchain4 + nchain5; + + // number of hits required per track + int minhits= 16; + + // first part of filenames + char filestub[] = "/raid4/nsj/sim/bggen/dedx_tree_g4_50932_bggen_50932"; + char filestub2[] = "/raid4/nsj/sim/bggen/mcr_nt_4120_tree_g4_50932_bggen"; + char filestub3[] = "/raid4/nsj/sim/dedx_nt_tree_g4_50932_bggen_1000000"; + char filestub4[] = "/raid4/nsj/sim/dedx_nt_tree_g4_50392_bggen_1000000"; + char filestub5[] = "/raid4/nsj/sim/dedx_nt_tree_g4_50932_bggen_1000000"; + + // set up the chain of the root files which will be processed in turn + + TChain chain("dedx_tree"); + + for (int i=0; i78) continue; + if (r>1) continue; + + // find which histo to fill + + int ihisto = 0; + + while (theta > thetamax[ihisto]) ihisto++; + + // cout << "theta " << theta << " ihisto " << ihisto << " thetamax[ihisto] " << thetamax[ihisto] << endl; + + + // fill the histograms + + if (ihisto < nhistos) { + if (charge > 0) { + dedxp_pos[ihisto]->Fill(p,dedx); + } else { + dedxp_neg[ihisto]->Fill(p,dedx); + } + } + + } + + + /* + + // Draw the histograms on the screen - VERY slow + + for (int i=30; i<31; i++) { //probably don't want 0 to nhistos unless nhistos is very small + + new TCanvas; + dedxp_pos[i]->Draw("colz"); + + new TCanvas; + dedxp_neg[i]->Draw("colz"); + + } + */ + + // create an output file & save the histograms into it + + TFile *f = new TFile(Form("sim_binned_histos_%ihits_%ifiles.root",minhits,totalnchain),"RECREATE"); + + for (int i=0; iWrite(); + for (int i=0; iWrite(); + // f->ls(); // list the file contents + f->Close(); + + +} diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/dots.C b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/dots.C new file mode 100644 index 0000000000..1797c70fcf --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/dots.C @@ -0,0 +1,291 @@ +{ + TTree* t = (TTree*)gDirectory->Get("t"); + + int nreal; + float realtheta[2000], realdedx[2000], realf[2000], realp[2000]; + int realpid[2000]; + + int nfake; + float faketheta[2000], fakededx[2000], fakef[2000], fakep[2000]; + int fakepid[2000]; + + int nall; + float alltheta[2000], alldedx[2000], allf[2000], allp[2000]; + int allpid[2000]; + + + + int np90; + float p90theta[2000], p90dedx[2000], p90f[2000], p90p[2000]; + + int np80; + float p80theta[2000], p80dedx[2000], p80f[2000], p80p[2000]; + + int np70; + float p70theta[2000], p70dedx[2000], p70f[2000], p70p[2000]; + + int np60; + float p60theta[2000], p60dedx[2000], p60f[2000], p60p[2000]; + + int np50; + float p50theta[2000], p50dedx[2000], p50f[2000], p50p[2000]; + + int np45; + float p45theta[2000], p45dedx[2000], p45f[2000], p45p[2000]; + + + int npi90; + float pi90theta[2000], pi90dedx[2000], pi90f[2000], pi90p[2000]; + + int npi80; + float pi80theta[2000], pi80dedx[2000], pi80f[2000], pi80p[2000]; + + int npi70; + float pi70theta[2000], pi70dedx[2000], pi70f[2000], pi70p[2000]; + + int npi60; + float pi60theta[2000], pi60dedx[2000], pi60f[2000], pi60p[2000]; + + int npi50; + float pi50theta[2000], pi50dedx[2000], pi50f[2000], pi50p[2000]; + + int npi45; + float pi45theta[2000], pi45dedx[2000], pi45f[2000], pi45p[2000]; + + + + + + float theta, dedx,f, p; + int pid; + + t->SetBranchAddress("f",&f); + t->SetBranchAddress("p",&p); + t->SetBranchAddress("pid",&pid); + t->SetBranchAddress("theta",&theta); + t->SetBranchAddress("dedx",&dedx); + + double n = t->GetEntries(); + + + nall = (int)n; + nreal = 0; + nfake = 0; + + for (int i=0; i< n; i++) { + t->GetEntry(i); + + alltheta[i] = theta; + alldedx[i] = dedx; + allpid[i] = pid; + allf[i] = f; + allp[i] = p; + + if (pid>0) { + realtheta[nreal] = theta; + realdedx[nreal] = dedx; + realpid[nreal] = pid; + realf[nreal] = f; + realp[nreal] = p; + nreal++; + } + + if (pid<0) { + faketheta[nfake] = theta; + fakededx[nfake] = dedx; + fakepid[nfake] = pid; + fakef[nfake] = f; + fakep[nfake] = p; + nfake++; + } + + if (pid<0) continue; + + if (pid == 1) { + + if (theta==90.0) { + p90theta[np90] = theta; + p90dedx[np90] = dedx; + p90f[np90] = f; + p90p[np90] = p; + np90++; + } else if (theta == 80.0) { + p80theta[np80] = theta; + p80dedx[np80] = dedx; + p80f[np80] = f; + p80p[np80] = p; + np80++; + } else if (theta == 70.0) { + p70theta[np70] = theta; + p70dedx[np70] = dedx; + p70f[np70] = f; + p70p[np70] = p; + np70++; + } else if (theta == 60.0) { + p60theta[np60] = theta; + p60dedx[np60] = dedx; + p60f[np60] = f; + p60p[np60] = p; + np60++; + } else if (theta == 50.0) { + p50theta[np50] = theta; + p50dedx[np50] = dedx; + p50f[np50] = f; + p50p[np50] = p; + np50++; + } else if (theta == 45.0) { + p45theta[np45] = theta; + p45dedx[np45] = dedx; + p45f[np45] = f; + p45p[np45] = p; + np45++; + } + + } else if (pid == 3) { + if (theta==90.0) { + pi90theta[npi90] = theta; + pi90dedx[npi90] = dedx; + pi90f[npi90] = f; + pi90p[npi90] = p; + npi90++; + } else if (theta == 80.0) { + pi80theta[npi80] = theta; + pi80dedx[npi80] = dedx; + pi80f[npi80] = f; + pi80p[npi80] = p; + npi80++; + } else if (theta == 70.0) { + pi70theta[npi70] = theta; + pi70dedx[npi70] = dedx; + pi70f[npi70] = f; + pi70p[npi70] = p; + npi70++; + } else if (theta == 60.0) { + pi60theta[npi60] = theta; + pi60dedx[npi60] = dedx; + pi60f[npi60] = f; + pi60p[npi60] = p; + npi60++; + } else if (theta == 50.0) { + pi50theta[npi50] = theta; + pi50dedx[npi50] = dedx; + pi50f[npi50] = f; + pi50p[npi50] = p; + npi50++; + } else if (theta == 45.0) { + pi45theta[npi45] = theta; + pi45dedx[npi45] = dedx; + pi45f[npi45] = f; + pi45p[npi45] = p; + npi45++; + } + } // end if pid is 1 or 3 + } + + gStyle->SetMarkerStyle(20); + gStyle->SetMarkerSize(0.5); + + TGraph *greal = new TGraph(nreal,realtheta,realdedx); + greal->SetMarkerColor(kMagenta); + + TGraph *gfake = new TGraph(nfake,faketheta,fakededx); + gfake->SetMarkerColor(38); + + TMultiGraph *mg = new TMultiGraph(); + mg->Add(greal); + mg->Add(gfake); + mg->SetTitle("Real (pink) and interpolated (grey) correction factors; theta (degrees); dE/dx (keV/cm)"); + mg->Draw("AP"); + + new TCanvas; + TGraph2D *gg = new TGraph2D(nall,alltheta,alldedx,allf); + gg->SetTitle("Correction factors; theta (degrees); dE/dx (keV/cm); correction factor"); + gg->Draw("APcolz"); + + return; + + int col[] = {900, 880, 860, 433, 819, 801, 805, 14}; + + TGraph *gp90 = new TGraph(np90,p90dedx,p90f); + TGraph *gp80 = new TGraph(np80,p80dedx,p80f); + TGraph *gp70 = new TGraph(np70,p70dedx,p70f); + TGraph *gp60 = new TGraph(np60,p60dedx,p60f); + TGraph *gp50 = new TGraph(np50,p50dedx,p50f); + TGraph *gp45 = new TGraph(np45,p45dedx,p45f); + + TGraph *gpi90 = new TGraph(npi90,pi90dedx,pi90f); + TGraph *gpi80 = new TGraph(npi80,pi80dedx,pi80f); + TGraph *gpi70 = new TGraph(npi70,pi70dedx,pi70f); + + gStyle->SetMarkerSize(0.7); + + gp90->SetMarkerSize(0.7); + gp80->SetMarkerSize(0.7); + gp70->SetMarkerSize(0.7); + gp60->SetMarkerSize(0.7); + gp50->SetMarkerSize(0.7); + gp45->SetMarkerSize(0.7); + + gpi90->SetMarkerSize(0.7); + gpi80->SetMarkerSize(0.7); + gpi70->SetMarkerSize(0.7); + + + gp90->SetMarkerColor(col[0]); + gp80->SetMarkerColor(col[1]); + gp70->SetMarkerColor(col[2]); + gp60->SetMarkerColor(col[3]); + gp50->SetMarkerColor(col[4]); + gp45->SetMarkerColor(col[5]); + + gpi90->SetMarkerColor(col[0]); + gpi80->SetMarkerColor(col[1]); + gpi70->SetMarkerColor(col[2]); + + + + TLegend *leg1 = new TLegend(0.92,0.3,0.99,0.9); + leg1->AddEntry(gp90,"90","P"); + leg1->AddEntry(gp80,"80","P"); + leg1->AddEntry(gp70,"70","P"); + leg1->AddEntry(gp60,"60","P"); + leg1->AddEntry(gp50,"50","P"); + leg1->AddEntry(gp45,"45","P"); + leg1->SetBorderSize(0); + + TLegend *leg2 = new TLegend(0.92,0.3,0.99,0.9); + leg2->AddEntry(gpi90,"90","P"); + leg2->AddEntry(gpi80,"80","P"); + leg2->AddEntry(gpi70,"70","P"); + leg2->SetBorderSize(0); + + + + + + new TCanvas; + + TMultiGraph *mgp = new TMultiGraph(); + mgp->Add(gp90); + mgp->Add(gp80); + mgp->Add(gp70); + mgp->Add(gp60); + mgp->Add(gp50); + mgp->Add(gp45); + mgp->SetTitle("protons;dE/dx;correction factor"); + mgp->Draw("AP"); + leg1->Draw(); + + new TCanvas; + + TMultiGraph *mgpi = new TMultiGraph(); + mgpi->Add(gpi90); + mgpi->Add(gpi80); + mgpi->Add(gpi70); + mgpi->SetTitle("pions;dE/dx;correction factor"); + + mgpi->Draw("AP"); + leg2->Draw(); + + +} diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/drawalldedxcuts.C b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/drawalldedxcuts.C new file mode 100644 index 0000000000..1b0664a5c7 --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/drawalldedxcuts.C @@ -0,0 +1,140 @@ + +// draws de/dx cut functions on top of the histogram specified from the file specified + +// draws these for every histogram and saves to folder cutpics + +{ + gStyle->SetMarkerStyle(20); + gStyle->SetMarkerSize(0.5); + gStyle->SetOptFit(1); + gStyle->SetOptStat(0); + gStyle->SetOptLogz(1); + + + char histofile[300]; + + sprintf(histofile, "/home/njarvis/online/dedx/binned_histos_20hits_150files.root"); + + TFile *f = new TFile(histofile,"READ"); + + if (!f) printf("Cannot find %s\n",histofile); + if (!f) gApplication->Terminate(); + + // functional forms for cut functions + + TF1 *fp = new TF1("fp","exp([0]*(x + [1])) + [2]", 0.0, 4.0); // between protons and deuterons + TF1 *fk = new TF1("fk","exp([0]*(x + [1])) + [2]", 0.0, 0.85); // between protons and kaons + TF1 *fpi = new TF1("fpi","exp([0]*(x + [1])) + [2]", 0.0, 0.425); // below pions + + + TCanvas *c = new TCanvas("c","c",1200,1200); + c->Divide(2,2); + + + TH2D* hpos; // histo pointer + + int ihisto; // which histo to draw (draws 2 per loop) + + + for (int i=80; i>10; i--) { + + ihisto = i; + + + + // plots are symmetric about 90 degrees (histo 80) + + // cout << "Difference between selected angle and 90 degrees is " << fabs(ihisto-80) << endl; + + + if (fabs(ihisto-80) < 10) { // ihisto 71-89 + fp->SetParameters(-4.7, -0.88, 4.0); // good for 90 degr + fk->SetParameters(-4.8,-0.65,2.3); // good for 90 degr + fpi->SetParameters(-6.5, -0.4, 2.4); + + } else if (fabs(ihisto-80) < 20) { // ihisto 61-70, 90-99 + fp->SetParameters(-4.6, -0.90, 4.8); // 80 degr + fk->SetParameters(-5.0,-0.64,3.0); + fpi->SetParameters(-6.5, -0.4, 2.8); + + } else if (fabs(ihisto-80) < 30) { // ihisto 51-60, 100-109 + fp->SetParameters(-4.6, -0.95, 4.8); // 70 degr + fk->SetParameters(-4.8,-0.69,2.7); + fpi->SetParameters(-6.5, -0.43, 2.75); + + } else if (fabs(ihisto-80) < 40) { // ihisto 41-50, 110-119 + fp->SetParameters(-4.6, -0.98, 5.2); // 60 degr + fk->SetParameters(-4.6,-0.7,2.3); // good at 60 degr + fpi->SetParameters(-6.5, -0.43, 2.75); + + } else if (fabs(ihisto-80) < 50) { // ihisto 31-40, 120-129 + fp->SetParameters(-4.6, -0.98, 5.2); // 50 degr + fk->SetParameters(-4.6,-0.69,2.7); // good at 50 degr + fpi->SetParameters(-6.5, -0.43, 2.75); + + } else { //if (fabs(ihisto-80) < 55) { // ihisto < 31, > 129 + fp->SetParameters(-4.6, -0.98, 5.2); // 45 degr same as 60 + fk->SetParameters(-5.0, -0.7, 2.9); + fpi->SetParameters(-6.5, -0.43, 2.75); + } + + + + + hpos = (TH2D*)gDirectory->Get(Form("dedxp_q+[%i]",i)); + + hpos->GetYaxis()->SetRangeUser(0,15); + hpos->GetXaxis()->SetRangeUser(0,2); + + c->cd(1); + + hpos->Draw("colz"); + + fp->Draw("same"); + fk->Draw("same"); + fpi->Draw("same"); + + + c->cd(2); + + gPad->SetLogz(0); + hpos->Draw("colz"); + + fp->Draw("same"); + fk->Draw("same"); + fpi->Draw("same"); + + + ihisto = 160 - i; // will be 80 for 80 (90 degr) nvm + + hpos = (TH2D*)gDirectory->Get(Form("dedxp_q+[%i]",ihisto)); + + hpos->GetYaxis()->SetRangeUser(0,15); + hpos->GetXaxis()->SetRangeUser(0,2); + + c->cd(3); + + hpos->Draw("colz"); + + fp->Draw("same"); + fk->Draw("same"); + fpi->Draw("same"); + + + c->cd(4); + + gPad->SetLogz(0); + hpos->Draw("colz"); + + fp->Draw("same"); + fk->Draw("same"); + fpi->Draw("same"); + + c->SaveAs(Form("cutpics/cuts_%.idegr.png",10+i)); + + } + +} + + + diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/drawdedxcuts.C b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/drawdedxcuts.C new file mode 100644 index 0000000000..c1e6daf8d7 --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/drawdedxcuts.C @@ -0,0 +1,83 @@ + +// draws de/dx cut functions on top of the histogram specified from the file specified + +{ + gStyle->SetMarkerStyle(20); + gStyle->SetMarkerSize(0.5); + gStyle->SetOptFit(1); + gStyle->SetOptStat(0); + gStyle->SetOptLogz(1); + + + int ihisto = 30; // draw histo # 30 (for 40 degrees, theta=ihisto+10) + + char histofile[300]; + + sprintf(histofile, "/home/njarvis/online/dedx/binned_histos_20hits_150files.root"); + + TFile *f = new TFile(histofile,"READ"); + + if (!f) printf("Cannot find %s\n",histofile); + if (!f) gApplication->Terminate(); + + + // functional forms for cut functions + + TF1 *fp = new TF1("fp","exp([0]*(x + [1])) + [2]", 0.0, 4.0); // between protons and deuterons + TF1 *fk = new TF1("fk","exp([0]*(x + [1])) + [2]", 0.0, 0.85); // between protons and kaons + TF1 *fpi = new TF1("fpi","exp([0]*(x + [1])) + [2]", 0.0, 0.425); // below pions + + + // plots are symmetric about 90 degrees (histo 80) + + cout << "Difference between selected angle and 90 degrees is " << fabs(ihisto-80) << endl; + + if (fabs(ihisto-80) < 10) { // ihisto 71-89 + fp->SetParameters(-4.7, -0.88, 4.0); // good for 90 degr + fk->SetParameters(-4.8,-0.65,2.3); // good for 90 degr + fpi->SetParameters(-6.5, -0.4, 2.4); + + } else if (fabs(ihisto-80) < 20) { // ihisto 61-70, 90-99 + fp->SetParameters(-4.6, -0.90, 4.8); // 80 degr + fk->SetParameters(-5.0,-0.64,3.0); + fpi->SetParameters(-6.5, -0.4, 2.8); + + } else if (fabs(ihisto-80) < 30) { // ihisto 51-60, 100-109 + fp->SetParameters(-4.6, -0.95, 4.8); // 70 degr + fk->SetParameters(-4.8,-0.69,2.7); + fpi->SetParameters(-6.5, -0.43, 2.75); + + } else if (fabs(ihisto-80) < 40) { // ihisto 41-50, 110-119 + fp->SetParameters(-4.6, -0.98, 5.2); // 60 degr + fk->SetParameters(-4.6,-0.7,2.3); // good at 60 degr + fpi->SetParameters(-6.5, -0.43, 2.75); + + } else if (fabs(ihisto-80) < 50) { // ihisto 31-40, 120-129 + fp->SetParameters(-4.6, -0.98, 5.2); // 50 degr + fk->SetParameters(-4.6,-0.69,2.7); // good at 50 degr + fpi->SetParameters(-6.5, -0.43, 2.75); + + } else { //if (fabs(ihisto-80) < 55) { // ihisto < 31, > 129 + fp->SetParameters(-4.6, -0.98, 5.2); // 45 degr same as 60 + fk->SetParameters(-5.0, -0.7, 2.9); + fpi->SetParameters(-6.5, -0.43, 2.75); + } + + + TH2D* hpos = (TH2D*)gDirectory->Get(Form("dedxp_q+[%i]",ihisto)); + + hpos->GetYaxis()->SetRangeUser(0,15); + hpos->GetXaxis()->SetRangeUser(0,3); + + new TCanvas; + hpos->Draw("colz"); + + fp->Draw("same"); + fk->Draw("same"); + fpi->Draw("same"); + + +} + + + diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/drawdedxfromtchain.C b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/drawdedxfromtchain.C new file mode 100644 index 0000000000..01103a1562 --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/drawdedxfromtchain.C @@ -0,0 +1,161 @@ +{ + + gStyle->SetMarkerStyle(20); + gStyle->SetMarkerSize(0.5); + gStyle->SetOptFit(1); + gStyle->SetOptStat(0); + gStyle->SetOptLogz(1); + + + + int minhits= 16; + + const float dedxmax = 20; + // char filename[500]; + // sprintf(filename, "/raid12/gluex/rawdata2/Run050932/dedx_tree_050932_000_015.root"); + // sprintf(filename, "/raid12/gluex/rawdata2/Run050766/dedx_moremass_tree_050766_000_010.root"); + + // char filestub1[] = "/raid12/gluex/rawdata2/Run050932/dedx_tree2_pi_d_tree_050932"; + // char filestub1[] = "/raid12/gluex/rawdata2/Run050932/dedx_tree_nt_fixv4_tree_050932"; + // char filestub1[] = "/raid4/nsj/rawdata/Run050932/dedx_tree_nt_v7_s5a_tree_050932"; + // char filestub1[] = "/raid4/nsj/rawdata/Run050932/dedx_tree_418_nt_v3_tree_050932"; + + const bool saveplots=1; + /* + char filestub1[] = "/raid2/nsj/rawdata/Run050932/dedx_tree_4_18_1_v7_s5a2_int_tree_050932"; + char title[] = "corrected data (v7, s5a2, recon 4.18.1, truncated)"; + char filename_suffix[] = "corrected_4_18_1_v7_s5a2"; +*/ + + /* + char filestub1[] = "/raid2/nsj/rawdata/Run050932/dedx_tree_4_18_1_tree_050932"; + char title[] = "original data (recon 4.18.1, truncated)"; + char filename_suffix[] = "4_18_1"; + */ + /* + char filestub1[] = "/raid4/nsj/sim/bggen/mcr_4120_tree_g4_50932_bggen_50932"; + char title[] = "MC data (sim 4.12.0, truncated)"; + char filename_suffix[] = "sim_4_12_0"; + */ + + int nchain1=4; //number of files to add to the tchain + char filestub1[] = "/raid2/nsj/rawdata/Run050932/dedx_tree_4_20_newmatrix_notrunc_050932"; + char title[] = "corrected data, without hit truncation"; + char filename_suffix[] = "4_20_050932_corr_notrunc"; + + + // char filestub1[] = "/raid2/nsj/rawdata/Run050932/dedx_tree_4_20_newmatrix_050932"; + // char title[] = "corrected data, with hit truncation"; + // char filename_suffix[] = "4_20_050932_corr"; + + // char filestub1[] = "/raid2/nsj/rawdata/Run050932/dedx_tree_4_20_nocorr_050932"; + // char title[] = "original data, with hit truncation"; + // char filename_suffix[] = "4_20_050932_nocorr"; + + // char filestub1[] = "/raid2/nsj/rawdata/Run050932/dedx_tree_4_20_nocorr_notrunc_050932"; + // char title[] = "original data, without hit truncation"; + // char filename_suffix[] = "4_20_050932_nocorr_notrunc"; + + + + TChain chain("dedx_tree"); + + for (int i=0; iGet("dedx_tree"); + + chain.SetBranchAddress("eventnumber",&eventnumber); + chain.SetBranchAddress("ctracks",&ctracks); + chain.SetBranchAddress("ndedxhits",&nhits); + chain.SetBranchAddress("x",&x); + chain.SetBranchAddress("y",&y); + chain.SetBranchAddress("z",&z); + chain.SetBranchAddress("r",&r); + chain.SetBranchAddress("phi",&phi); + chain.SetBranchAddress("theta",&theta); + chain.SetBranchAddress("charge",&charge); + chain.SetBranchAddress("p",&p); + chain.SetBranchAddress("dedx",&dedx); + chain.SetBranchAddress("dedx_int",&dedx_int); + + int n = (int)chain.GetEntries(); + + if (n<1) cout << "empty tree\n"; + if (n<1) gApplication->Terminate(); + + cout << n << " entries\n"; + + + TH2D *dedxp = new TH2D("dedxp",Form("%s;p (GeV/c);dE/dx (keV/cm)",title),300,0,3,600,0,dedxmax); + TH2D *dedxtheta = new TH2D("dedxtheta",Form("%s;theta (degrees);dE/dx (keV/cm)",title),180,0,180,600,0,dedxmax); + + TH2D *intdedxp = new TH2D("intdedxp",Form("%s;p (GeV/c);dE/dx from integral (keV/cm)",title),300,0,3,600,0,dedxmax); + TH2D *intdedxtheta = new TH2D("intdedxtheta",Form("%s;theta (degrees);dE/dx from integral (keV/cm)",title),180,0,180,600,0,dedxmax); + + + for (int i=0; i78) continue; + if (r>1) continue; + + // find which histo to fill + + + dedxp->Fill(p,dedx); + dedxtheta->Fill(theta,dedx); + intdedxp->Fill(p,dedx_int); + intdedxtheta->Fill(theta,dedx_int); + + } + + if (saveplots) { + + /* + TFile *f = new TFile(Form("dedx_%s.root",filename_suffix),"RECREATE"); + + f->cd(); + + dedxp->Write(); + dedxtheta->Write(); + intdedxp->Write(); + intdedxtheta->Write(); + */ + new TCanvas; + gStyle->SetOptStat(0); + dedxp->Draw("colz"); + gPad->SetLogz(); + // dedxp->Draw("colz"); + gPad->SaveAs(Form("dedxp_%s.png",filename_suffix)); + new TCanvas; + dedxtheta->Draw("colz"); + gPad->SaveAs(Form("dedxtheta_%s.png",filename_suffix)); + /* new TCanvas; + intdedxp->Draw("colz"); + gPad->SaveAs(Form("intdedxp_%s.png",filename_suffix)); + new TCanvas; + intdedxtheta->Draw("colz"); + gPad->SaveAs(Form("intdedxtheta_%s.png",filename_suffix)); */ + } +} diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/drawsurf2.C b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/drawsurf2.C new file mode 100644 index 0000000000..6f41d209e8 --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/drawsurf2.C @@ -0,0 +1,138 @@ +{ + + const float dpmax = 0; //0.5; // exclude d with mtm > dpmax + + TTree* t = (TTree*)gDirectory->Get("t"); + + int nentries = (int)t->GetEntries(); + float theta, dedx, f, p; + int pid; + t->SetBranchAddress("dedx",&dedx); + t->SetBranchAddress("theta",&theta); + t->SetBranchAddress("f",&f); + t->SetBranchAddress("pid",&pid); + t->SetBranchAddress("p",&p); + + + // theta 21 to 137 + // dedx 2.30382 to 15.0024 + + TCanvas *c = new TCanvas("c","Graph2D example",0,0,800,600); + Double_t x, y, z; + // Int_t np = 200; + TGraph2D *gtot = new TGraph2D(); + gtot->SetNpx(116); + gtot->SetNpy(254); + + TGraph2D *gprotons = new TGraph2D(); + gprotons->SetNpx(116); + gprotons->SetNpy(254); + + TGraph2D *gpions = new TGraph2D(); + gpions->SetNpx(116); + gpions->SetNpy(254); + + TGraph2D *gdeuterons = new TGraph2D(); + gdeuterons->SetNpx(116); + gdeuterons->SetNpy(254); + + TGraph2D *gcombo = new TGraph2D(); + gcombo->SetNpx(116); + gcombo->SetNpy(254); + + TGraph2D *gfakes = new TGraph2D(); + gfakes->SetNpx(116); + gfakes->SetNpy(254); + + int np=0; + int nf=0; + int nd=0; + int npi=0; + int nc=0; + + int nt=0; // use for drawing more than one particle, just draw over one in a diff colour. + // to get the axis range right + + for (Int_t N=0; NGetEntry(N); + + if (dedx>10) continue; + + x = theta; //2*P*(r->Rndm(N))-P; + y = dedx; //2*P*(r->Rndm(N))-P; + z = f; //(sin(x)/x)*(sin(y)/y)+0.2; + if (pid==-1) gfakes->SetPoint(nf,x,y,z); + if (pid==0 && p < dpmax) gdeuterons->SetPoint(nd,x,y,z); + if (pid==1) gprotons->SetPoint(np,x,y,z); + if (pid==3) gpions->SetPoint(npi,x,y,z); + if (pid==4) gcombo->SetPoint(nc,x,y,z); + + if (fabs(pid)!=0 && pid!=4) gtot->SetPoint(nt,x,y,z); + if (fabs(pid) != 0 && pid !=4) nt++; + + if (pid==0 && pSetPoint(nt,x,y,z); + if (pid==0 && pSetPalette(1); + gtot->SetTitle("dE/dx correction factor;theta;dedx;f"); + gtot->SetMarkerColor(0); + gtot->Draw("surf1"); + + new TCanvas; + + // gtot->GetYaxis()->SetRangeUser(2,10); + + gtot->Draw("surf1"); + + gprotons->SetMarkerStyle(20); + gpions->SetMarkerStyle(20); + gcombo->SetMarkerStyle(20); + gdeuterons->SetMarkerStyle(20); + gfakes->SetMarkerStyle(20); + + gprotons->SetMarkerSize(0.4); + gprotons->SetMarkerColor(kRed); + gprotons->Draw("psame"); + + gpions->SetMarkerColor(kBlue); + gpions->SetMarkerSize(0.4); + gpions->Draw("psame"); + + gcombo->SetMarkerColor(kMagenta); + gcombo->SetMarkerSize(0.4); + // gcombo->Draw("psame"); + + gdeuterons->SetMarkerColor(906); + gdeuterons->SetMarkerSize(0.5); + gdeuterons->Draw("psame"); + + gfakes->SetMarkerColor(38); + gfakes->SetMarkerSize(0.5); + gfakes->Draw("psame"); + + // g + + // return c; + + new TCanvas; + gtot->Draw("p"); + gprotons->Draw("psame"); + gpions->Draw("psame"); + gcombo->Draw("psame"); + gdeuterons->Draw("psame"); + gfakes->Draw("psame"); + /* + TH2D *h = dt->GetHistogram(); + + new TCanvas; + h->Draw("surf"); + */ +} diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/make_histo_into_text_file.C b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/make_histo_into_text_file.C new file mode 100644 index 0000000000..fbef1dd450 --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/make_histo_into_text_file.C @@ -0,0 +1,210 @@ +{ + + TTree* t = (TTree*)gDirectory->Get("t"); + + int nentries = (int)t->GetEntries(); + float theta, dedx, f, p; + int pid; + t->SetBranchAddress("dedx",&dedx); + t->SetBranchAddress("theta",&theta); + t->SetBranchAddress("f",&f); + t->SetBranchAddress("pid",&pid); + t->SetBranchAddress("p",&p); + + const int THETA_MIN = 21; + const int THETA_MAX = 136; + int ntheta = 1 + THETA_MAX - THETA_MIN ; + + // theta 22 to 135, extended to 21 to 136 + // dedx 2.30382 to 15.0024, extended to 20 + // histo dedx 2.30740 to 20 + + + TCanvas *c = new TCanvas("c","Graph2D example",0,0,800,600); + Double_t x, y, z; + // Int_t np = 200; + + const int npy = 1+ 100; + + TGraph2D *gtot = new TGraph2D(); + gtot->SetNpx(ntheta); + gtot->SetNpy(npy); + + TGraph2D *gprotons = new TGraph2D(); + gprotons->SetNpx(ntheta); + gprotons->SetNpy(npy); + + TGraph2D *gpions = new TGraph2D(); + gpions->SetNpx(ntheta); + gpions->SetNpy(npy); + + TGraph2D *gdeuterons = new TGraph2D(); + gdeuterons->SetNpx(ntheta); + gdeuterons->SetNpy(npy); + + TGraph2D *gcombo = new TGraph2D(); + gcombo->SetNpx(ntheta); + gcombo->SetNpy(npy); + + TGraph2D *gfakes = new TGraph2D(); + gfakes->SetNpx(ntheta); + gfakes->SetNpy(npy); + + int np=0; + int nf=0; + int nd=0; + int npi=0; + int nc=0; + + int nt=0; // use for drawing more than one particle, just draw over one in a diff colour. + // to get the axis range right + + float dcut = 0.46; + + for (Int_t N=0; NGetEntry(N); + x = theta; //2*P*(r->Rndm(N))-P; + y = dedx; //2*P*(r->Rndm(N))-P; + z = f; //(sin(x)/x)*(sin(y)/y)+0.2; + if (pid==-1) gfakes->SetPoint(nf,x,y,z); + if (pid==0 && p < dcut) gdeuterons->SetPoint(nd,x,y,z); + if (pid==1) gprotons->SetPoint(np,x,y,z); + if (pid==3) gpions->SetPoint(npi,x,y,z); + if (pid==4) gcombo->SetPoint(nc,x,y,z); + + if (fabs(pid)!=0 && pid!=4) gtot->SetPoint(nt,x,y,z); + if (fabs(pid) != 0 && pid !=4) nt++; + + // if (pid==0 && pSetPoint(nt,x,y,z); + // if (pid==0 && pSetPalette(1); + gtot->SetTitle("dE/dx correction factor;theta;dedx;f"); + gtot->SetMarkerColor(0); + gtot->Draw("surf1"); + + TH2D* histo = gtot->GetHistogram(); + + new TCanvas; + histo->Draw("surf1"); + + cout << npy << endl; + + char filename[100]; + sprintf(filename,"surface_histo_%i_%i.png",ntheta,npy); + + c1->SaveAs(filename); + + sprintf(filename,"surface_histo_%i_%i.root",ntheta,npy); + + TFile* newfile = new TFile(filename,"RECREATE"); + histo->Write(); + newfile->Close(); + + // histo has npy bins in y and ntheta bins in x + + // bins are ordered underflow data overflow + // ie 2 extra bins per histo + + double thetamin = histo->GetXaxis()->GetXmin(); + double thetamax = histo->GetXaxis()->GetXmax(); + + double dedxmin = histo->GetYaxis()->GetXmin(); + double dedxmax = histo->GetYaxis()->GetXmax(); + + double thetastep = (thetamax - thetamin) / (float)(ntheta - 1); + double dedxstep = (dedxmax - dedxmin) / (float)(npy - 1); + + sprintf(filename,"dedx_corrections_%i_%i.txt",ntheta,npy); + FILE *outfile = fopen(filename,"w"); + fprintf(outfile,"%i values of theta\n",ntheta); //histo->GetNbinsX()); + fprintf(outfile,"%f min theta\n",thetamin); + fprintf(outfile,"%f max theta\n",thetamax); + fprintf(outfile,"%f theta step\n",thetastep); + + fprintf(outfile,"%i values of dedx\n",npy); + fprintf(outfile,"%f min dedx\n",dedxmin); + fprintf(outfile,"%f max dedx\n",dedxmax); + fprintf(outfile,"%f dedx step\n",dedxstep); + + fprintf(outfile,"\n"); + + + int totalbins = (ntheta+2)*(npy+2); + int ix = 0; + + for (int ibin = ntheta+2; ibin0 && ix<=ntheta) fprintf(outfile,"%f\n",histo->GetBinContent(ibin)); + ix++; + if (ix== ntheta+2) ix=0; + } + + fclose(outfile); + + /* + int ibin=0; + + cout << histo->GetNbinsX() << "x bins" << endl; + cout << histo->GetNbinsY() << "y bins" << endl; + + cout << endl; + + cout << "first bin 0\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + ibin=ntheta-1; + cout << "last bin in bottom row " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + + + ibin = ntheta; + cout << "bin ntheta " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + ibin = ntheta+1; + cout << "bin ntheta+1: " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + ibin = ntheta+2; + cout << "bin ntheta+2: " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + + ibin = ntheta+3; + cout << "bin ntheta+3: " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + + ibin = ntheta+4; + cout << "bin ntheta+4: " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + + // cout << histo->GetXaxis()->GetBinLowEdge(ibin); + // cout << histo->GetYaxis()->GetBinLowEdge(ibin); + // cout << histo->GetBinContent(ibin); + + + ibin = (ntheta+0)*npy; + cout << "last bin " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + + ibin = (ntheta+1)*npy; + cout << "last overflow bin " << ibin << "\n"; + cout << histo->GetXaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetYaxis()->GetBinLowEdge(ibin) << "\t" << histo->GetBinContent(ibin) << "\n"; + + + +*/ + +} diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/plotdata.C b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/plotdata.C new file mode 100644 index 0000000000..e0e1e48bee --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/plotdata.C @@ -0,0 +1,177 @@ +{ + + // plots peak positions from root tree in fitresults.root or fittedsimresults.root + + gStyle->SetMarkerStyle(20); + gStyle->SetMarkerSize(0.4); + + const int fixrange=0; //draw invisible graph to set the axis range + + int plotpid=1; // 1=protons, 3=pions + + const int sim=1; + + char treefile[50]; + sprintf(treefile, "fitresults.root"); + if (sim) sprintf(treefile, "fittedsimresults.root"); + + TFile *f = new TFile(treefile,"READ"); + if (f) cout << "Opened " << treefile << endl; + if (!f) exit; + + TTree *t = (TTree*)gDirectory->Get("t"); + + + char pidname[5][30] = {"deuteron","proton","K+","pi+","proton and pi+"}; + + int pid; + float p, theta; + double nentries = t->GetEntries(); + + t->SetBranchAddress("pid",&pid); + t->SetBranchAddress("p",&p); + t->SetBranchAddress("theta",&theta); + + + // t->SetBranchAddress("n",nentries); + + // use arrays so that all pids can be included more easily + + // 0: d 1: p 2: K 3: pi + + double ht, mean, sigma; + double err_c, err_m, err_s; + + t->SetBranchAddress("ht",&ht); + t->SetBranchAddress("mean",&mean); + t->SetBranchAddress("sigma",&sigma); + t->SetBranchAddress("err_h",&err_c); + t->SetBranchAddress("err_m",&err_m); + t->SetBranchAddress("err_s",&err_s); + + const int np = 55; // number of momentum values + const int nt = 161; // number of theta values + + + int n[np]={0}; + double d[np][nt], th[np][nt]; + double pmin; + double pset[np]; + int ngraph; + + + + double pstart=0; + double pstop=5.0; + + if (plotpid==1) pstart = 0.3; + + // if (plotpid==3) pstart = 0.15; + if (plotpid==3) pstop=0.46; // the fits from 0.5 + are not good + + double minht; + + if (plotpid==1) minht=30; //30 + if (plotpid==1 && sim==1) minht = 25; + + if (plotpid==3) minht=100; + + for (int i=0; i<(int)nentries; i++) { + + t->GetEntry(i); + + if (pid != plotpid) continue; + if (p < pstart) continue; + if (p > pstop) continue; + if (ht <= minht) continue; + + //printf("theta %.0f p %.1f p_mean %.1f \n",theta,p,mean[1]); + + ngraph = -1; + + for (int j=0; j= 0) break; + + } + + + int npoints = n[ngraph]; + + d[ngraph][npoints] = mean; + th[ngraph][npoints] = theta; + + n[ngraph]++; + + } + + + int col[] = {1, 2, 609, 880, 860, 66, 834, 426, 800, 95 , 804, 892, 852, 872, 38}; + + TGraph *g[np]; + + TMultiGraph *mg = new TMultiGraph(); + + TGraph *gghost; + if (fixrange) { + double ymin = 0, ymax = 15.5; + if (plotpid==3) ymin = 2; + if (plotpid==3) ymax = 4; + double xghost[] = {45,130}; + double yghost[] = {ymin,ymax}; + gghost = new TGraph(2,xghost,yghost); + gghost->SetMarkerColor(0); + mg->Add(gghost); + } + + + + TLegend *leg = new TLegend(0.91,0.1, 0.999, 0.9); + leg->SetBorderSize(0); + leg->SetMargin(0); + TLegendEntry *le[55]; + + le[0] = leg->AddEntry(g[0],"GeV/c",""); + le[0]->SetTextColor(1); + + + for (int i=0; i<15; i++) { + + + if (n[i]>0) { + + g[i] = new TGraph(n[i],th[i],d[i]); + g[i]->SetMarkerColor(col[i]); + mg->Add(g[i]); + float thisp = pstart + i*0.05; + le[i+1] = leg->AddEntry(g[i],Form("%.2f",thisp),""); + le[i+1]->SetTextColor(col[i]); + + } + + } + + char title[300] = ""; + char simlabel[8]; + char peaktype[8]; + + if (sim) sprintf(simlabel,"bggen "); + + if (plotpid==1) sprintf(peaktype,"mean"); + if (plotpid==3) sprintf(peaktype,"mode"); + + sprintf(title,"%s of %s%s dE/dx;theta (degrees);dE/dx (keV/cm)", peaktype, simlabel, pidname[plotpid]); + + mg->SetTitle(title); + + + mg->Draw("AP"); + + // if (plotpid==1) mg->GetYaxis()->SetRangeUser(0,15.5); + // if (plotpid==3) mg->GetYaxis()->SetRangeUser(2,3.6); + + leg->Draw(); + +} diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/plotf.C b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/plotf.C new file mode 100644 index 0000000000..61a0eb3d27 --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/plotf.C @@ -0,0 +1,169 @@ +{ + + // mtm range is 0.3 to 3 + + gStyle->SetMarkerStyle(20); + gStyle->SetMarkerSize(0.4); + + + const int fixrange=0; //draw invisible graph to set the axis range + + int plotpid = 1; + + char treefile[50]; + + sprintf(treefile, "dedx_correction_factors.root"); + + + TFile *f = new TFile(treefile,"READ"); + if (f) cout << "Opened " << treefile << endl; + if (!f) exit; + + TTree *t = (TTree*)gDirectory->Get("t"); + + + char pidname[5][30] = {"deuteron","proton","K+","pi+","proton and pi+"}; + + int pid; + float p, theta, dedx, cf; + double nentries = t->GetEntries(); + + t->SetBranchAddress("pid",&pid); + t->SetBranchAddress("p",&p); + t->SetBranchAddress("f",&cf); + t->SetBranchAddress("theta",&theta); + t->SetBranchAddress("dedx",&dedx); + + const int np = 15; // number of momentum values + const int nt = 86; // number of theta values + + + int n[np]={0}; + double d[np][nt], th[np][nt]; + double pmin; + double pset[np]; + int ngraph; + + + + double pstart=0; + double pstop=1.0; + + if (plotpid==1) pstart = 0.35; + + if (plotpid==3) pstart = 0.15; + if (plotpid==3) pstop=0.46; // the fits from 0.5 + are not good + + + for (int i=0; i<(int)nentries; i++) { + + t->GetEntry(i); + + if (pid != plotpid) continue; + // if (p < pstart) continue; + // if (p > pstop) continue; + + //printf("theta %.0f p %.1f p_mean %.1f \n",theta,p,mean[1]); + + ngraph = -1; + + for (int j=0; j= 0) break; + + } + + //cout << " p " << p << " graph # " << ngraph << endl; + + int npoints = n[ngraph]; + + d[ngraph][npoints] = cf; //dedx; + th[ngraph][npoints] = theta; + + + //cout << "ngraph " << n[ngraph] << " dedx " << d[ngraph][n[ngraph]] << " theta " << th[ngraph][n[ngraph]]<< endl; + + + n[ngraph]++; + + } + + + int col[] = {1, 2, 609, 880, 860, 66, 834, 426, 800, 95 , 804, 892, 852, 872, 38}; + + TGraph *g[np]; + + TMultiGraph *mg = new TMultiGraph(); + + TGraph *gghost; + if (fixrange) { + double ymin = 0, ymax = 15.5; + if (plotpid==3) ymin = 2; + if (plotpid==3) ymax = 4; + double xghost[] = {45,130}; + double yghost[] = {ymin,ymax}; + gghost = new TGraph(2,xghost,yghost); + gghost->SetMarkerColor(0); + mg->Add(gghost); + } + + + + TLegend *leg = new TLegend(0.91,0.1, 0.999, 0.9); + leg->SetBorderSize(0); + leg->SetMargin(0); + TLegendEntry *le[55]; + + le[0] = leg->AddEntry(g[0],"GeV/c",""); + le[0]->SetTextColor(1); + + + for (int i=0; i<15; i++) { + + + if (n[i]>0) { + + g[i] = new TGraph(n[i],th[i],d[i]); + g[i]->SetMarkerColor(col[i]); + mg->Add(g[i]); + float thisp = pstart + i*0.05; + le[i+1] = leg->AddEntry(g[i],Form("%.2f",thisp),""); + le[i+1]->SetTextColor(col[i]); + + } + + } + + char title[300] = ""; + + sprintf(title,"%s dE/dx correction factors;theta (degrees);correction factor", pidname[plotpid]); + + mg->SetTitle(title); + +// mg->SetTitle(Form("mean of %s dE/dx band;theta (degrees);dE/dx (keV/cm)", pidname[plotpid])); + + //if (ntracks>0) mg->SetTitle(Form("%i track events, mean of %s dE/dx band;theta (degrees);dE/dx (keV/cm)",ntracks, pidname[plotpid])); + + mg->Draw("AP"); + + // if (plotpid==1) mg->GetYaxis()->SetRangeUser(0,15.5); + // if (plotpid==3) mg->GetYaxis()->SetRangeUser(2,3.6); + + leg->Draw(); + + + + +/* + gStyle->SetMarkerSize(0.2); + gStyle->SetMarkerStyle(20); + + t->Draw("p_mean:theta>>h(150,0,3,1700,0,17)","!p_fitstat && p<0.31"); + + h->Draw(); + */ + + +} diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/plotfdedx.C b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/plotfdedx.C new file mode 100644 index 0000000000..7cfb282089 --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/plottingscripts/plotfdedx.C @@ -0,0 +1,158 @@ +{ + + // mtm range is 0.3 to 3 + + gStyle->SetMarkerStyle(20); + gStyle->SetMarkerSize(0.4); + + + const int fixrange=0; //draw invisible graph to set the axis range + + int plotpid = 1; + + char treefile[50]; + + sprintf(treefile, "dedx_correction_factors.root"); + + + TFile *f = new TFile(treefile,"READ"); + if (f) cout << "Opened " << treefile << endl; + if (!f) exit; + + TTree *t = (TTree*)gDirectory->Get("t"); + + + char pidname[5][30] = {"deuteron","proton","K+","pi+","proton and pi+"}; + + int pid; + float p, theta, dedx, cf; + double nentries = t->GetEntries(); + + t->SetBranchAddress("pid",&pid); + t->SetBranchAddress("p",&p); + t->SetBranchAddress("f",&cf); + t->SetBranchAddress("theta",&theta); + t->SetBranchAddress("dedx",&dedx); + + const int maxnp = 713; // + const int maxnpi = 357; // = 2*15; // make this npid//15; // number of momentum values + const int nt = 86; // number of theta values + + + + int np = 0 ; + int npi = 0; + double pd[maxnp], pth[maxnp]; + double piond[maxnpi], pionth[maxnpi]; + + double pmin; + + int ngraph; + + + + double pstart=0; + double pstop=1.0; + + if (plotpid==1) pstart = 0.35; + + if (plotpid==3) pstart = 0.15; + if (plotpid==3) pstop=0.46; // the fits from 0.5 + are not good + + + for (int i=0; i<(int)nentries; i++) { + + t->GetEntry(i); + + // if (pid != plotpid) continue; + + pstart = 0.34; + pstop = 0.86; + + if (pid==3) pstart = 0.14; + if (pid==3) pstop=0.46; // the fits from 0.5 + are not good + + if (p < pstart) continue; + if (p > pstop) continue; + + //printf("theta %.0f p %.1f p_mean %.1f \n",theta,p,mean[1]); + + ngraph = 0; + if (pid==3) ngraph = 1; + + //cout << " p " << p << " graph # " << ngraph << endl; + + if (pid==1) { + pd[np] = dedx; + pth[np] = theta; + np++; + } else if (pid == 3) { + piond[npi] = dedx; + pionth[npi] = theta; + npi++; + } + + + } + + + int col[] = {1, 2, 609, 880, 860, 66, 834, 426, 800, 95 , 804, 892, 852, 872, 38}; + + + + TMultiGraph *mg = new TMultiGraph(); + + TGraph *gghost; + if (fixrange) { + double ymin = 0, ymax = 15.5; + if (plotpid==3) ymin = 2; + if (plotpid==3) ymax = 4; + double xghost[] = {45,130}; + double yghost[] = {ymin,ymax}; + gghost = new TGraph(2,xghost,yghost); + gghost->SetMarkerColor(0); + mg->Add(gghost); + } + + + TGraph* gp = new TGraph(np,pth,pd); + gp->SetMarkerColor(col[0]); + mg->Add(gp); + + TGraph* gpi = new TGraph(npi,pionth,piond); + gpi->SetMarkerColor(col[1]); + mg->Add(gpi); + + + TLegend *leg = new TLegend(0.91,0.1, 0.999, 0.9); + leg->SetBorderSize(0); + leg->SetMargin(0); + TLegendEntry *le[55]; + + le[0] = leg->AddEntry(gp,"protons",""); + le[0]->SetTextColor(col[0]); + + le[1] = leg->AddEntry(gpi,"pions",""); + le[1]->SetTextColor(col[1]); + + char title[300] = ""; + + sprintf(pidname[plotpid],""); + sprintf(title,"%s dE/dx values used for correction factors;theta (degrees);dE/dx (keV/cm)", pidname[plotpid]); + + mg->SetTitle(title); + +// mg->SetTitle(Form("mean of %s dE/dx band;theta (degrees);dE/dx (keV/cm)", pidname[plotpid])); + + //if (ntracks>0) mg->SetTitle(Form("%i track events, mean of %s dE/dx band;theta (degrees);dE/dx (keV/cm)",ntracks, pidname[plotpid])); + + mg->Draw("AP"); + + // if (plotpid==1) mg->GetYaxis()->SetRangeUser(0,15.5); + // if (plotpid==3) mg->GetYaxis()->SetRangeUser(2,3.6); + + leg->Draw(); + + + +} diff --git a/src/plugins/Utilities/retired/dedx_tree/scripts/smooth_histo.C b/src/plugins/Utilities/retired/dedx_tree/scripts/smooth_histo.C new file mode 100644 index 0000000000..5279d85ee4 --- /dev/null +++ b/src/plugins/Utilities/retired/dedx_tree/scripts/smooth_histo.C @@ -0,0 +1,201 @@ +{ + + TTree* t = (TTree*)gDirectory->Get("t"); + + int nentries = (int)t->GetEntries(); + float theta, dedx, f, p; + int pid; + t->SetBranchAddress("dedx",&dedx); + t->SetBranchAddress("theta",&theta); + t->SetBranchAddress("f",&f); + t->SetBranchAddress("pid",&pid); + t->SetBranchAddress("p",&p); + + const int THETA_MIN = 21; + const int THETA_MAX = 136; + int ntheta = 1 + THETA_MAX - THETA_MIN ; + + // theta 22 to 135, extended to 21 to 136 + // dedx 2.30382 to 15.0024, extended to 20 + // histo dedx 2.30740 to 20 + + + TCanvas *c = new TCanvas("c","Graph2D example",0,0,800,600); + Double_t x, y, z; + // Int_t np = 200; + + const int npy = 1+ 100; + + TGraph2D *gtot = new TGraph2D(); + gtot->SetNpx(ntheta); + gtot->SetNpy(npy); + + TGraph2D *gprotons = new TGraph2D(); + gprotons->SetNpx(ntheta); + gprotons->SetNpy(npy); + + TGraph2D *gpions = new TGraph2D(); + gpions->SetNpx(ntheta); + gpions->SetNpy(npy); + + TGraph2D *gdeuterons = new TGraph2D(); + gdeuterons->SetNpx(ntheta); + gdeuterons->SetNpy(npy); + + TGraph2D *gcombo = new TGraph2D(); + gcombo->SetNpx(ntheta); + gcombo->SetNpy(npy); + + TGraph2D *gfakes = new TGraph2D(); + gfakes->SetNpx(ntheta); + gfakes->SetNpy(npy); + + int np=0; + int nf=0; + int nd=0; + int npi=0; + int nc=0; + + int nt=0; // use for drawing more than one particle, just draw over one in a diff colour. + // to get the axis range right + + float dcut = 0.46; + + for (Int_t N=0; NGetEntry(N); + x = theta; //2*P*(r->Rndm(N))-P; + y = dedx; //2*P*(r->Rndm(N))-P; + z = f; //(sin(x)/x)*(sin(y)/y)+0.2; + if (pid==-1) gfakes->SetPoint(nf,x,y,z); + if (pid==0 && p < dcut) gdeuterons->SetPoint(nd,x,y,z); + if (pid==1) gprotons->SetPoint(np,x,y,z); + if (pid==3) gpions->SetPoint(npi,x,y,z); + if (pid==4) gcombo->SetPoint(nc,x,y,z); + + if (fabs(pid)!=0 && pid!=4) gtot->SetPoint(nt,x,y,z); + if (fabs(pid) != 0 && pid !=4) nt++; + + // if (pid==0 && pSetPoint(nt,x,y,z); + // if (pid==0 && pSetPalette(1); + gtot->SetTitle("dE/dx correction factor;theta;dedx;f"); + gtot->SetMarkerColor(0); + gtot->Draw("surf1"); + + TH2D* histo = gtot->GetHistogram(); + + new TCanvas; + histo->Draw("surf1"); + gPad->SaveAs("v6_original.png"); + + /* + TH2D *histo3 = (TH2D*)histo->Clone(); + histo3->Smooth(1,"3a"); + + new TCanvas; + histo3->Draw("surf1"); + gPad->SaveAs("v6_smoothed3a.png"); + + TH2D *histo5 = (TH2D*)histo->Clone(); + histo5->Smooth(1,"5a"); + + new TCanvas; + histo5->Draw("surf1"); + gPad->SaveAs("v6_smoothed5a.png"); + + TH2D *histo5b = (TH2D*)histo->Clone(); + histo5b->Smooth(1,"5b"); + + new TCanvas; + histo5b->Draw("surf1"); + gPad->SaveAs("v6_smoothed5b.png"); + */ + + //cout << npy << endl; + + char smoothopt[3]; + + sprintf(smoothopt,"5a"); + + histo->Smooth(1,smoothopt); + histo->Smooth(1,smoothopt); + + sprintf(smoothopt,"5a_twice"); + + char filename[100]; + sprintf(filename,"surface_histo_s%s_%i_%i.png",smoothopt,ntheta,npy); + + gPad->SaveAs(filename); + + sprintf(filename,"surface_histo_s%s_%i_%i.root",smoothopt,ntheta,npy); + + + // + + + + + // + TFile* newfile = new TFile(filename,"RECREATE"); + histo->Write(); + newfile->Close(); + + new TCanvas; + histo->Draw("colz"); + gPad->SaveAs(Form("surface_histo_flat_s%s_%i_%i.png",smoothopt,ntheta,npy)); + + + // histo has npy bins in y and ntheta bins in x + + // bins are ordered underflow data overflow + // ie 2 extra bins per histo + + double thetamin = histo->GetXaxis()->GetXmin(); + double thetamax = histo->GetXaxis()->GetXmax(); + + double dedxmin = histo->GetYaxis()->GetXmin(); + double dedxmax = histo->GetYaxis()->GetXmax(); + + double thetastep = (thetamax - thetamin) / (float)(ntheta - 1); + double dedxstep = (dedxmax - dedxmin) / (float)(npy - 1); + + sprintf(filename,"dedx_corrections_s%s_%i_%i.txt",smoothopt,ntheta,npy); + + cout << "Writing text file " << filename << endl; + + FILE *outfile = fopen(filename,"w"); + fprintf(outfile,"%i values of theta\n",ntheta); //histo->GetNbinsX()); + fprintf(outfile,"%f min theta\n",thetamin); + fprintf(outfile,"%f max theta\n",thetamax); + fprintf(outfile,"%f theta step\n",thetastep); + + fprintf(outfile,"%i values of dedx\n",npy); + fprintf(outfile,"%f min dedx\n",dedxmin); + fprintf(outfile,"%f max dedx\n",dedxmax); + fprintf(outfile,"%f dedx step\n",dedxstep); + + fprintf(outfile,"\n"); + + + int totalbins = (ntheta+2)*(npy+2); + int ix = 0; + + for (int ibin = ntheta+2; ibin0 && ix<=ntheta) fprintf(outfile,"%f\n",histo->GetBinContent(ibin)); + ix++; + if (ix== ntheta+2) ix=0; + } + + fclose(outfile); + + +} diff --git a/src/plugins/Utilities/dumpcandidates/JEventProcessor_dumpcandidates.cc b/src/plugins/Utilities/retired/dumpcandidates/JEventProcessor_dumpcandidates.cc similarity index 100% rename from src/plugins/Utilities/dumpcandidates/JEventProcessor_dumpcandidates.cc rename to src/plugins/Utilities/retired/dumpcandidates/JEventProcessor_dumpcandidates.cc diff --git a/src/plugins/Utilities/dumpcandidates/JEventProcessor_dumpcandidates.h b/src/plugins/Utilities/retired/dumpcandidates/JEventProcessor_dumpcandidates.h similarity index 100% rename from src/plugins/Utilities/dumpcandidates/JEventProcessor_dumpcandidates.h rename to src/plugins/Utilities/retired/dumpcandidates/JEventProcessor_dumpcandidates.h diff --git a/src/plugins/Utilities/dumpthrowns/SConscript b/src/plugins/Utilities/retired/dumpcandidates/SConscript similarity index 100% rename from src/plugins/Utilities/dumpthrowns/SConscript rename to src/plugins/Utilities/retired/dumpcandidates/SConscript diff --git a/src/plugins/Utilities/dumpthrowns/JEventProcessor_dumpthrowns.cc b/src/plugins/Utilities/retired/dumpthrowns/JEventProcessor_dumpthrowns.cc similarity index 100% rename from src/plugins/Utilities/dumpthrowns/JEventProcessor_dumpthrowns.cc rename to src/plugins/Utilities/retired/dumpthrowns/JEventProcessor_dumpthrowns.cc diff --git a/src/plugins/Utilities/dumpthrowns/JEventProcessor_dumpthrowns.h b/src/plugins/Utilities/retired/dumpthrowns/JEventProcessor_dumpthrowns.h similarity index 100% rename from src/plugins/Utilities/dumpthrowns/JEventProcessor_dumpthrowns.h rename to src/plugins/Utilities/retired/dumpthrowns/JEventProcessor_dumpthrowns.h diff --git a/src/plugins/Utilities/l3bdt/SConscript b/src/plugins/Utilities/retired/dumpthrowns/SConscript similarity index 100% rename from src/plugins/Utilities/l3bdt/SConscript rename to src/plugins/Utilities/retired/dumpthrowns/SConscript diff --git a/src/plugins/Utilities/es_test/JEventProcessor_es_test.cc b/src/plugins/Utilities/retired/es_test/JEventProcessor_es_test.cc similarity index 100% rename from src/plugins/Utilities/es_test/JEventProcessor_es_test.cc rename to src/plugins/Utilities/retired/es_test/JEventProcessor_es_test.cc diff --git a/src/plugins/Utilities/es_test/JEventProcessor_es_test.h b/src/plugins/Utilities/retired/es_test/JEventProcessor_es_test.h similarity index 100% rename from src/plugins/Utilities/es_test/JEventProcessor_es_test.h rename to src/plugins/Utilities/retired/es_test/JEventProcessor_es_test.h diff --git a/src/plugins/Utilities/es_test/SConscript b/src/plugins/Utilities/retired/es_test/SConscript similarity index 100% rename from src/plugins/Utilities/es_test/SConscript rename to src/plugins/Utilities/retired/es_test/SConscript diff --git a/src/plugins/Utilities/l3bdt/JEventProcessor_L3BDTtree.cc b/src/plugins/Utilities/retired/l3bdt/JEventProcessor_L3BDTtree.cc similarity index 100% rename from src/plugins/Utilities/l3bdt/JEventProcessor_L3BDTtree.cc rename to src/plugins/Utilities/retired/l3bdt/JEventProcessor_L3BDTtree.cc diff --git a/src/plugins/Utilities/l3bdt/JEventProcessor_L3BDTtree.h b/src/plugins/Utilities/retired/l3bdt/JEventProcessor_L3BDTtree.h similarity index 100% rename from src/plugins/Utilities/l3bdt/JEventProcessor_L3BDTtree.h rename to src/plugins/Utilities/retired/l3bdt/JEventProcessor_L3BDTtree.h diff --git a/src/plugins/Utilities/l3bdt/README b/src/plugins/Utilities/retired/l3bdt/README similarity index 100% rename from src/plugins/Utilities/l3bdt/README rename to src/plugins/Utilities/retired/l3bdt/README diff --git a/src/plugins/Utilities/retired/l3bdt/SConscript b/src/plugins/Utilities/retired/l3bdt/SConscript new file mode 100644 index 0000000000..3005e0e91d --- /dev/null +++ b/src/plugins/Utilities/retired/l3bdt/SConscript @@ -0,0 +1,13 @@ + + +import sbms + +# get env object and clone it +Import('*') +env = env.Clone() + +sbms.AddDANA(env) +sbms.AddROOT(env) +sbms.plugin(env) + + diff --git a/src/plugins/Utilities/l3bdt/trainBDT.C b/src/plugins/Utilities/retired/l3bdt/trainBDT.C similarity index 100% rename from src/plugins/Utilities/l3bdt/trainBDT.C rename to src/plugins/Utilities/retired/l3bdt/trainBDT.C diff --git a/src/plugins/Utilities/rawevent/JEventProcessor_rawevent.cc b/src/plugins/Utilities/retired/rawevent/JEventProcessor_rawevent.cc similarity index 100% rename from src/plugins/Utilities/rawevent/JEventProcessor_rawevent.cc rename to src/plugins/Utilities/retired/rawevent/JEventProcessor_rawevent.cc diff --git a/src/plugins/Utilities/rawevent/JEventProcessor_rawevent.h b/src/plugins/Utilities/retired/rawevent/JEventProcessor_rawevent.h similarity index 100% rename from src/plugins/Utilities/rawevent/JEventProcessor_rawevent.h rename to src/plugins/Utilities/retired/rawevent/JEventProcessor_rawevent.h diff --git a/src/plugins/Utilities/rawevent/Makefile b/src/plugins/Utilities/retired/rawevent/Makefile similarity index 100% rename from src/plugins/Utilities/rawevent/Makefile rename to src/plugins/Utilities/retired/rawevent/Makefile diff --git a/src/plugins/Utilities/rawevent/SConscript b/src/plugins/Utilities/retired/rawevent/SConscript similarity index 100% rename from src/plugins/Utilities/rawevent/SConscript rename to src/plugins/Utilities/retired/rawevent/SConscript diff --git a/src/plugins/Utilities/rawevent/makeTranslationTable.py b/src/plugins/Utilities/retired/rawevent/makeTranslationTable.py similarity index 100% rename from src/plugins/Utilities/rawevent/makeTranslationTable.py rename to src/plugins/Utilities/retired/rawevent/makeTranslationTable.py diff --git a/src/plugins/Utilities/rawevent/mc2coda.c b/src/plugins/Utilities/retired/rawevent/mc2coda.c similarity index 100% rename from src/plugins/Utilities/rawevent/mc2coda.c rename to src/plugins/Utilities/retired/rawevent/mc2coda.c diff --git a/src/plugins/Utilities/rawevent/mc2coda.h b/src/plugins/Utilities/retired/rawevent/mc2coda.h similarity index 100% rename from src/plugins/Utilities/rawevent/mc2coda.h rename to src/plugins/Utilities/retired/rawevent/mc2coda.h diff --git a/src/plugins/Utilities/rawevent/mc2coda_modules.h b/src/plugins/Utilities/retired/rawevent/mc2coda_modules.h similarity index 100% rename from src/plugins/Utilities/rawevent/mc2coda_modules.h rename to src/plugins/Utilities/retired/rawevent/mc2coda_modules.h diff --git a/src/plugins/Utilities/rawevent/mc2coda_random.cc b/src/plugins/Utilities/retired/rawevent/mc2coda_random.cc similarity index 100% rename from src/plugins/Utilities/rawevent/mc2coda_random.cc rename to src/plugins/Utilities/retired/rawevent/mc2coda_random.cc diff --git a/src/plugins/monitoring/SConscript b/src/plugins/monitoring/SConscript index abec710145..025fc0852f 100644 --- a/src/plugins/monitoring/SConscript +++ b/src/plugins/monitoring/SConscript @@ -5,32 +5,24 @@ subdirs = [ 'occupancy_online', 'lowlevel_online', 'BCAL_online', -'BCAL_LEDonline', -'BCAL_LED', -'BCAL_LED_time', 'CCAL_online', 'CDC_dedx', 'CDC_online', +'CDC_roc_hits', 'DAQ_online', 'fa125_itrig', 'FCAL_online', 'FCAL_invmass', -'FCAL_cpp', 'FDC_online', 'FDC_Efficiency', 'FMWPC_online', 'FMWPC_Performance', 'PSC_online', 'RF_online', -'pedestal_online', 'PS_online', 'PSPair_online', 'PS_flux', 'ST_online_lowlevel', -'ST_online_tracking', -'ST_online_Tresolution', -'ST_online_efficiency', -'ST_online_multi', 'TAGGER_online', 'TAGH_online', 'TAGM_online', @@ -38,42 +30,28 @@ subdirs = [ 'TOF_online', 'TRIG_online', 'TPOL_online', -'CDC_expert', 'CDC_expert_2', 'TOF_TDC_shift', 'BCAL_Eff', 'BCAL_inv_mass', 'CDC_drift', -'CDC_PerStrawReco', -'CDC_roc_hits', -'EPICS_dump', 'CDC_Efficiency', 'L1_online', 'highlevel_online', 'BCAL_Hadronic_Eff', 'FCAL_Hadronic_Eff', -'SC_Eff','TOF_Eff', -'TS_scaler', 'TrackingPulls', -'TrackingPulls_straight', -'timing_online', 'DIRC_online', -'TRD_online', -'RSAI_KO', 'BEAM_online', 'lumi_mon', -'scaler_primex' -,'cppFMWPC' -,'cppFMWPC_ana' -,'cpp_itrig' -,'primex-online' -,'HELI_online' +'HELI_online' ] -#'L3_online', -#'CODA_online', -#'EVNT_online', +optdirs = ['TAGH_doubles','ST_ZEff', 'CODA_online', 'EVNT_online'] +optdirs.extend(['BCAL_LED_time', 'BCAL_LED', 'FCAL_cpp', 'pedestal_online', 'CDC_PerStrawReco' ]) +optdirs.extend(['EPICS_dump', 'SC_Eff','TOF_Eff', 'TS_scaler', 'TrackingPulls_straight', 'timing_online', 'TRD_online' ]) +optdirs.extend(['RSAI_KO', 'scaler_primex', 'cppFMWPC', 'cppFMWPC_ana', 'cpp_itrig', 'primex-online' ]) +sbms.OptionallyBuild(env, optdirs) -sbms.OptionallyBuild(env, ['TAGH_doubles','ST_ZEff']) SConscript(dirs=subdirs, exports='env osname', duplicate=0) diff --git a/src/plugins/monitoring/BCAL_LED_time/JEventProcessor_BCAL_LED_time.cc b/src/plugins/monitoring/retired/BCAL_LED_time/JEventProcessor_BCAL_LED_time.cc similarity index 100% rename from src/plugins/monitoring/BCAL_LED_time/JEventProcessor_BCAL_LED_time.cc rename to src/plugins/monitoring/retired/BCAL_LED_time/JEventProcessor_BCAL_LED_time.cc diff --git a/src/plugins/monitoring/BCAL_LED_time/JEventProcessor_BCAL_LED_time.h b/src/plugins/monitoring/retired/BCAL_LED_time/JEventProcessor_BCAL_LED_time.h similarity index 100% rename from src/plugins/monitoring/BCAL_LED_time/JEventProcessor_BCAL_LED_time.h rename to src/plugins/monitoring/retired/BCAL_LED_time/JEventProcessor_BCAL_LED_time.h diff --git a/src/plugins/monitoring/BCAL_LED_time/SConscript b/src/plugins/monitoring/retired/BCAL_LED_time/SConscript similarity index 100% rename from src/plugins/monitoring/BCAL_LED_time/SConscript rename to src/plugins/monitoring/retired/BCAL_LED_time/SConscript diff --git a/src/plugins/monitoring/BCAL_LED_time/scripts/hist_read.C b/src/plugins/monitoring/retired/BCAL_LED_time/scripts/hist_read.C similarity index 100% rename from src/plugins/monitoring/BCAL_LED_time/scripts/hist_read.C rename to src/plugins/monitoring/retired/BCAL_LED_time/scripts/hist_read.C diff --git a/src/plugins/monitoring/BCAL_LED_time/scripts/hist_read.py b/src/plugins/monitoring/retired/BCAL_LED_time/scripts/hist_read.py similarity index 100% rename from src/plugins/monitoring/BCAL_LED_time/scripts/hist_read.py rename to src/plugins/monitoring/retired/BCAL_LED_time/scripts/hist_read.py diff --git a/src/plugins/monitoring/BCAL_LED_time/scripts/plot_diff_const.py b/src/plugins/monitoring/retired/BCAL_LED_time/scripts/plot_diff_const.py similarity index 100% rename from src/plugins/monitoring/BCAL_LED_time/scripts/plot_diff_const.py rename to src/plugins/monitoring/retired/BCAL_LED_time/scripts/plot_diff_const.py diff --git a/src/plugins/monitoring/CDC_expert/JEventProcessor_CDC_expert.cc b/src/plugins/monitoring/retired/CDC_expert/JEventProcessor_CDC_expert.cc similarity index 100% rename from src/plugins/monitoring/CDC_expert/JEventProcessor_CDC_expert.cc rename to src/plugins/monitoring/retired/CDC_expert/JEventProcessor_CDC_expert.cc diff --git a/src/plugins/monitoring/CDC_expert/JEventProcessor_CDC_expert.h b/src/plugins/monitoring/retired/CDC_expert/JEventProcessor_CDC_expert.h similarity index 100% rename from src/plugins/monitoring/CDC_expert/JEventProcessor_CDC_expert.h rename to src/plugins/monitoring/retired/CDC_expert/JEventProcessor_CDC_expert.h diff --git a/src/plugins/monitoring/CDC_expert/SConscript b/src/plugins/monitoring/retired/CDC_expert/SConscript similarity index 100% rename from src/plugins/monitoring/CDC_expert/SConscript rename to src/plugins/monitoring/retired/CDC_expert/SConscript diff --git a/src/plugins/monitoring/ST_online_Tresolution/JEventProcessor_ST_online_Tresolution.cc b/src/plugins/monitoring/retired/ST_online_Tresolution/JEventProcessor_ST_online_Tresolution.cc similarity index 100% rename from src/plugins/monitoring/ST_online_Tresolution/JEventProcessor_ST_online_Tresolution.cc rename to src/plugins/monitoring/retired/ST_online_Tresolution/JEventProcessor_ST_online_Tresolution.cc diff --git a/src/plugins/monitoring/ST_online_Tresolution/JEventProcessor_ST_online_Tresolution.h b/src/plugins/monitoring/retired/ST_online_Tresolution/JEventProcessor_ST_online_Tresolution.h similarity index 100% rename from src/plugins/monitoring/ST_online_Tresolution/JEventProcessor_ST_online_Tresolution.h rename to src/plugins/monitoring/retired/ST_online_Tresolution/JEventProcessor_ST_online_Tresolution.h diff --git a/src/plugins/monitoring/ST_online_Tresolution/SConscript b/src/plugins/monitoring/retired/ST_online_Tresolution/SConscript similarity index 100% rename from src/plugins/monitoring/ST_online_Tresolution/SConscript rename to src/plugins/monitoring/retired/ST_online_Tresolution/SConscript diff --git a/src/plugins/monitoring/ST_online_Tresolution/st_time_resolution.C b/src/plugins/monitoring/retired/ST_online_Tresolution/st_time_resolution.C similarity index 100% rename from src/plugins/monitoring/ST_online_Tresolution/st_time_resolution.C rename to src/plugins/monitoring/retired/ST_online_Tresolution/st_time_resolution.C diff --git a/src/plugins/monitoring/ST_online_efficiency/JEventProcessor_ST_online_efficiency.cc b/src/plugins/monitoring/retired/ST_online_efficiency/JEventProcessor_ST_online_efficiency.cc similarity index 100% rename from src/plugins/monitoring/ST_online_efficiency/JEventProcessor_ST_online_efficiency.cc rename to src/plugins/monitoring/retired/ST_online_efficiency/JEventProcessor_ST_online_efficiency.cc diff --git a/src/plugins/monitoring/ST_online_efficiency/JEventProcessor_ST_online_efficiency.h b/src/plugins/monitoring/retired/ST_online_efficiency/JEventProcessor_ST_online_efficiency.h similarity index 100% rename from src/plugins/monitoring/ST_online_efficiency/JEventProcessor_ST_online_efficiency.h rename to src/plugins/monitoring/retired/ST_online_efficiency/JEventProcessor_ST_online_efficiency.h diff --git a/src/plugins/monitoring/ST_online_efficiency/SConscript b/src/plugins/monitoring/retired/ST_online_efficiency/SConscript similarity index 100% rename from src/plugins/monitoring/ST_online_efficiency/SConscript rename to src/plugins/monitoring/retired/ST_online_efficiency/SConscript diff --git a/src/plugins/monitoring/ST_online_efficiency/ST_Monitoring_Eff.C b/src/plugins/monitoring/retired/ST_online_efficiency/ST_Monitoring_Eff.C similarity index 100% rename from src/plugins/monitoring/ST_online_efficiency/ST_Monitoring_Eff.C rename to src/plugins/monitoring/retired/ST_online_efficiency/ST_Monitoring_Eff.C diff --git a/src/plugins/monitoring/ST_online_multi/JEventProcessor_ST_online_multi.cc b/src/plugins/monitoring/retired/ST_online_multi/JEventProcessor_ST_online_multi.cc similarity index 100% rename from src/plugins/monitoring/ST_online_multi/JEventProcessor_ST_online_multi.cc rename to src/plugins/monitoring/retired/ST_online_multi/JEventProcessor_ST_online_multi.cc diff --git a/src/plugins/monitoring/ST_online_multi/JEventProcessor_ST_online_multi.h b/src/plugins/monitoring/retired/ST_online_multi/JEventProcessor_ST_online_multi.h similarity index 100% rename from src/plugins/monitoring/ST_online_multi/JEventProcessor_ST_online_multi.h rename to src/plugins/monitoring/retired/ST_online_multi/JEventProcessor_ST_online_multi.h diff --git a/src/plugins/monitoring/ST_online_multi/SConscript b/src/plugins/monitoring/retired/ST_online_multi/SConscript similarity index 100% rename from src/plugins/monitoring/ST_online_multi/SConscript rename to src/plugins/monitoring/retired/ST_online_multi/SConscript diff --git a/src/plugins/monitoring/ST_online_multi/st_multi.C b/src/plugins/monitoring/retired/ST_online_multi/st_multi.C similarity index 100% rename from src/plugins/monitoring/ST_online_multi/st_multi.C rename to src/plugins/monitoring/retired/ST_online_multi/st_multi.C diff --git a/src/plugins/monitoring/ST_online_tracking/JEventProcessor_ST_online_tracking.cc b/src/plugins/monitoring/retired/ST_online_tracking/JEventProcessor_ST_online_tracking.cc similarity index 100% rename from src/plugins/monitoring/ST_online_tracking/JEventProcessor_ST_online_tracking.cc rename to src/plugins/monitoring/retired/ST_online_tracking/JEventProcessor_ST_online_tracking.cc diff --git a/src/plugins/monitoring/ST_online_tracking/JEventProcessor_ST_online_tracking.h b/src/plugins/monitoring/retired/ST_online_tracking/JEventProcessor_ST_online_tracking.h similarity index 100% rename from src/plugins/monitoring/ST_online_tracking/JEventProcessor_ST_online_tracking.h rename to src/plugins/monitoring/retired/ST_online_tracking/JEventProcessor_ST_online_tracking.h diff --git a/src/plugins/monitoring/ST_online_tracking/SConscript b/src/plugins/monitoring/retired/ST_online_tracking/SConscript similarity index 100% rename from src/plugins/monitoring/ST_online_tracking/SConscript rename to src/plugins/monitoring/retired/ST_online_tracking/SConscript diff --git a/src/plugins/monitoring/ST_online_tracking/ST_Monitoring_2D.C b/src/plugins/monitoring/retired/ST_online_tracking/ST_Monitoring_2D.C similarity index 100% rename from src/plugins/monitoring/ST_online_tracking/ST_Monitoring_2D.C rename to src/plugins/monitoring/retired/ST_online_tracking/ST_Monitoring_2D.C diff --git a/src/plugins/monitoring/ST_online_tracking/ST_Monitoring_Eff.C b/src/plugins/monitoring/retired/ST_online_tracking/ST_Monitoring_Eff.C similarity index 100% rename from src/plugins/monitoring/ST_online_tracking/ST_Monitoring_Eff.C rename to src/plugins/monitoring/retired/ST_online_tracking/ST_Monitoring_Eff.C diff --git a/src/plugins/monitoring/ST_online_tracking/ST_Monitoring_Pid.C b/src/plugins/monitoring/retired/ST_online_tracking/ST_Monitoring_Pid.C similarity index 100% rename from src/plugins/monitoring/ST_online_tracking/ST_Monitoring_Pid.C rename to src/plugins/monitoring/retired/ST_online_tracking/ST_Monitoring_Pid.C