/* * evaluateResponsesNew.C * * Created on: Jun 11, 2020 * Author: John Russell * * Similar to "evaluateResponses.C", but for use with the OvR Fisher responses. */ #include "AnitaDataset.h" #include "AnitaTMVA.h" void evaluateResponsesNew(const char * sampleSumFileName, const char * sampleKLCoeffsFileName, const char * sampleResponseFileName) { // Open sample summary file and set branch addresses. TFile sampleSumFile(sampleSumFileName); TTree * sampleSumTree = (TTree *) sampleSumFile.Get("sampleA4"); sampleSumTree -> GetEntry(0); // Necessary to pass references in the following two lines. AnitaEventSummary * sampleSum = 0; sampleSumTree -> SetBranchAddress("summary", & sampleSum); int & run = sampleSum -> run; unsigned int & eventNumber = sampleSum -> eventNumber; int numEntries = sampleSumTree -> GetEntries(); // Open KL coefficient file and set branch addresses. TFile sampleKLCoeffsFile(sampleKLCoeffsFileName); // Structure from which coefficients stored, and related quantites. struct { double KLCoeffs[2][5]; // First index corresponds to if polarity flipped, second index to principal componennt order. double KLMag[2]; // Magnitude over KLCoeffs, index correspond to polarity flipped. unsigned int KLCoeff0Idx; // Based on sise of KLCoeffs[][0], 1 when KLCoeffs[1][0] > KLCoeffs[0][0], 0 otherwise. unsigned int KLMagIdx; // Based on size of KLMag, 1 when KLMag[1] > KLMag[0], 0 otherwise. Should correspond to polarity flip when KLMagIdx == 1. } coh, deconv, both; TTree * sampleKLCoeffsTree = (TTree *) sampleKLCoeffsFile.Get("KLCoeffsTree"); sampleKLCoeffsTree -> SetBranchAddress("coh", & coh); sampleKLCoeffsTree -> SetBranchAddress("deconv", & deconv); sampleKLCoeffsTree -> SetBranchAddress("both", & both); // Create TFile in which to place TTree produced by this macro. TFile sampleResponseFile(sampleResponseFileName, "recreate"); // Create new TTree to be added to the file, containing the evaluated responses. struct { double WAISHPol; double WAISVPol; double HiCal2A; double HiCal2B; double iceMC; double RCP; double LCP; double payloadBlast; double sun; double thermal; } Peng, Andrew, KLCoeffs, PengCDF, AndrewCDF, KLCoeffsCDF; // Create tree to save response values. TString responseBranchString = "WAISHPol/D:WAISVPol:HiCal2A:HiCal2B:iceMC:LCP:RCP:payloadBlast:sun:thermal"; TTree responseTree("responseTree", "TTree containing Fisher response values."); responseTree.Branch("run", & run); responseTree.Branch("eventNumber", & eventNumber); responseTree.Branch("Peng", & Peng, responseBranchString); responseTree.Branch("Andrew", & Andrew, responseBranchString); responseTree.Branch("KLCoeffs", & KLCoeffs, responseBranchString); responseTree.Branch("PengCDF", & PengCDF, responseBranchString); responseTree.Branch("AndrewCDF", & AndrewCDF, responseBranchString); responseTree.Branch("KLCoeffsCDF", & KLCoeffsCDF, responseBranchString); // Load in the data into appropriate TMVA readers. TMVA::Reader readerPeng, readerAndrew, readerKLCoeffs; // Add spectator variables. readerPeng.AddSpectator("run", & run); readerAndrew.AddSpectator("run", & run); readerKLCoeffs.AddSpectator("run", & run); int eventNumberInt; readerPeng.AddSpectator("eventNumber", & eventNumberInt); readerAndrew.AddSpectator("eventNumber", & eventNumberInt); readerKLCoeffs.AddSpectator("eventNumber", & eventNumberInt); int cohKLMagIdx, deconvKLMagIdx, bothKLMagIdx; readerKLCoeffs.AddSpectator("coh.KLMagIdx", & cohKLMagIdx); readerKLCoeffs.AddSpectator("deconv.KLMagIdx", & deconvKLMagIdx); readerKLCoeffs.AddSpectator("both.KLMagIdx", & bothKLMagIdx); // Add variables to be read in. float deconvImpulsivity; // Ordering of variables must reflect input file ordering. readerPeng.AddVariable("mostImpulsiveDeconvolvedFiltered(2).impulsivityMeasure", & deconvImpulsivity); readerAndrew.AddVariable("mostImpulsiveDeconvolvedFiltered(2).impulsivityMeasure", & deconvImpulsivity); float cohImpulsivity, deconvLinPolFrac, cohLinPolFrac, impulsivityDiff, peakHilbertDiff, deconvInstLinPolFrac, cohInstLinPolFrac, deconvPWG; readerAndrew.AddVariable("mostImpulsiveCoherentFiltered(2).impulsivityMeasure", & cohImpulsivity); readerAndrew.AddVariable("mostImpulsiveDeconvolvedFiltered(2).linearPolFrac()", & deconvLinPolFrac); readerAndrew.AddVariable("mostImpulsiveCoherentFiltered(2).linearPolFrac()", & cohLinPolFrac); readerAndrew.AddVariable("pow(mostImpulsiveDeconvolvedFiltered(2).impulsivityMeasure, 2) - pow(mostImpulsiveCoherentFiltered(2).impulsivityMeasure, 2)", & impulsivityDiff); readerAndrew.AddVariable("mostImpulsiveDeconvolvedFiltered(2).peakHilbert - mostImpulsiveCoherentFiltered(2).peakHilbert", & peakHilbertDiff); readerAndrew.AddVariable("mostImpulsiveDeconvolvedFiltered(2).instantaneousLinearPolFrac()", & deconvInstLinPolFrac); readerAndrew.AddVariable("mostImpulsiveCoherentFiltered(2).instantaneousLinearPolFrac()", & cohInstLinPolFrac); readerAndrew.AddVariable("mostImpulsiveDeconvolvedFiltered(2).fracPowerWindowGradient()", & deconvPWG); float cohKLMag, deconvKLMag, bothKLMag; readerKLCoeffs.AddVariable("coh.KLMag[coh.KLMagIdx]", & cohKLMag); readerKLCoeffs.AddVariable("deconv.KLMag[deconv.KLMagIdx]", & deconvKLMag); readerKLCoeffs.AddVariable("both.KLMag[both.KLMagIdx]", & bothKLMag); // Book MVA methods. Has to be done after adding variables. const char * resultsPath = "/home/nfs/jwruss/trainResponses/fisherOvRResponses/%sResults/weights/%sFisherOvRFactory_%s.weights.xml"; // Apparently ROOT doesn't know to equate "~" with "/home/nfs/jwruss" here? readerPeng.BookMVA("WAISHPol", TString::Format(resultsPath, "Peng", "Peng", "WAISHPol")); readerPeng.BookMVA("WAISVPol", TString::Format(resultsPath, "Peng", "Peng", "WAISVPol")); readerPeng.BookMVA("HiCal2A", TString::Format(resultsPath, "Peng", "Peng", "HiCal2A")); readerPeng.BookMVA("HiCal2A", TString::Format(resultsPath, "Peng", "Peng", "HiCal2B")); readerPeng.BookMVA("iceMC", TString::Format(resultsPath, "Peng", "Peng", "iceMC")); readerPeng.BookMVA("RCP", TString::Format(resultsPath, "Peng", "Peng", "RCP")); readerPeng.BookMVA("LCP", TString::Format(resultsPath, "Peng", "Peng", "LCP")); readerPeng.BookMVA("payloadBlast", TString::Format(resultsPath, "Peng", "Peng", "payloadBlast")); readerPeng.BookMVA("sun", TString::Format(resultsPath, "Peng", "Peng", "sun")); readerPeng.BookMVA("thermal", TString::Format(resultsPath, "Peng", "Peng", "thermal")); readerAndrew.BookMVA("WAISHPol", TString::Format(resultsPath, "Andrew", "Andrew", "WAISHPol")); readerAndrew.BookMVA("WAISVPol", TString::Format(resultsPath, "Andrew", "Andrew", "WAISVPol")); readerAndrew.BookMVA("HiCal2A", TString::Format(resultsPath, "Andrew", "Andrew", "HiCal2A")); readerAndrew.BookMVA("HiCal2A", TString::Format(resultsPath, "Andrew", "Andrew", "HiCal2B")); readerAndrew.BookMVA("iceMC", TString::Format(resultsPath, "Andrew", "Andrew", "iceMC")); readerAndrew.BookMVA("RCP", TString::Format(resultsPath, "Andrew", "Andrew", "RCP")); readerAndrew.BookMVA("LCP", TString::Format(resultsPath, "Andrew", "Andrew", "LCP")); readerAndrew.BookMVA("payloadBlast", TString::Format(resultsPath, "Andrew", "Andrew", "payloadBlast")); readerAndrew.BookMVA("sun", TString::Format(resultsPath, "Andrew", "Andrew", "sun")); readerAndrew.BookMVA("thermal", TString::Format(resultsPath, "Andrew", "Andrew", "thermal")); readerKLCoeffs.BookMVA("WAISHPol", TString::Format(resultsPath, "KLCoeffs", "KLCoeffs", "WAISHPol")); readerKLCoeffs.BookMVA("WAISVPol", TString::Format(resultsPath, "KLCoeffs", "KLCoeffs", "WAISVPol")); readerKLCoeffs.BookMVA("HiCal2A", TString::Format(resultsPath, "KLCoeffs", "KLCoeffs", "HiCal2A")); readerKLCoeffs.BookMVA("HiCal2A", TString::Format(resultsPath, "KLCoeffs", "KLCoeffs", "HiCal2B")); readerKLCoeffs.BookMVA("iceMC", TString::Format(resultsPath, "KLCoeffs", "KLCoeffs", "iceMC")); readerKLCoeffs.BookMVA("RCP", TString::Format(resultsPath, "KLCoeffs", "KLCoeffs", "RCP")); readerKLCoeffs.BookMVA("LCP", TString::Format(resultsPath, "KLCoeffs", "KLCoeffs", "LCP")); readerKLCoeffs.BookMVA("payloadBlast", TString::Format(resultsPath, "KLCoeffs", "KLCoeffs", "payloadBlast")); readerKLCoeffs.BookMVA("sun", TString::Format(resultsPath, "KLCoeffs", "KLCoeffs", "sun")); readerKLCoeffs.BookMVA("thermal", TString::Format(resultsPath, "KLCoeffs", "KLCoeffs", "thermal")); // Get TMVA::MethodFisher (insofar TMVA::IMethod relates to TMVA::MethodBase) in order to use rarity for signal type, as opposed to background type. TMVA::MethodFisher * PengWAISHPol = dynamic_cast(readerPeng.FindMVA("WAISHPol")); TMVA::MethodFisher * PengWAISVPol = dynamic_cast(readerPeng.FindMVA("WAISVPol")); TMVA::MethodFisher * PengHiCal2A = dynamic_cast(readerPeng.FindMVA("HiCal2A")); TMVA::MethodFisher * PengHiCal2B = dynamic_cast(readerPeng.FindMVA("HiCal2B")); TMVA::MethodFisher * PengIceMC = dynamic_cast(readerPeng.FindMVA("iceMC")); TMVA::MethodFisher * PengRCP = dynamic_cast(readerPeng.FindMVA("RCP")); TMVA::MethodFisher * PengLCP = dynamic_cast(readerPeng.FindMVA("LCP")); TMVA::MethodFisher * PengPayloadBlast = dynamic_cast(readerPeng.FindMVA("payloadBlast")); TMVA::MethodFisher * PengSun = dynamic_cast(readerPeng.FindMVA("sun")); TMVA::MethodFisher * PengThermal = dynamic_cast(readerPeng.FindMVA("thermal")); TMVA::MethodFisher * AndrewWAISHPol = dynamic_cast(readerAndrew.FindMVA("WAISHPol")); TMVA::MethodFisher * AndrewWAISVPol = dynamic_cast(readerAndrew.FindMVA("WAISVPol")); TMVA::MethodFisher * AndrewHiCal2A = dynamic_cast(readerAndrew.FindMVA("HiCal2A")); TMVA::MethodFisher * AndrewHiCal2B = dynamic_cast(readerAndrew.FindMVA("HiCal2B")); TMVA::MethodFisher * AndrewIceMC = dynamic_cast(readerAndrew.FindMVA("iceMC")); TMVA::MethodFisher * AndrewRCP = dynamic_cast(readerAndrew.FindMVA("RCP")); TMVA::MethodFisher * AndrewLCP = dynamic_cast(readerAndrew.FindMVA("LCP")); TMVA::MethodFisher * AndrewPayloadBlast = dynamic_cast(readerAndrew.FindMVA("payloadBlast")); TMVA::MethodFisher * AndrewSun = dynamic_cast(readerAndrew.FindMVA("sun")); TMVA::MethodFisher * AndrewThermal = dynamic_cast(readerAndrew.FindMVA("thermal")); TMVA::MethodFisher * KLCoeffsWAISHPol = dynamic_cast(readerKLCoeffs.FindMVA("WAISHPol")); TMVA::MethodFisher * KLCoeffsWAISVPol = dynamic_cast(readerKLCoeffs.FindMVA("WAISVPol")); TMVA::MethodFisher * KLCoeffsHiCal2A = dynamic_cast(readerKLCoeffs.FindMVA("HiCal2A")); TMVA::MethodFisher * KLCoeffsHiCal2B = dynamic_cast(readerKLCoeffs.FindMVA("HiCal2B")); TMVA::MethodFisher * KLCoeffsIceMC = dynamic_cast(readerKLCoeffs.FindMVA("iceMC")); TMVA::MethodFisher * KLCoeffsRCP = dynamic_cast(readerKLCoeffs.FindMVA("RCP")); TMVA::MethodFisher * KLCoeffsLCP = dynamic_cast(readerKLCoeffs.FindMVA("LCP")); TMVA::MethodFisher * KLCoeffsPayloadBlast = dynamic_cast(readerKLCoeffs.FindMVA("payloadBlast")); TMVA::MethodFisher * KLCoeffsSun = dynamic_cast(readerKLCoeffs.FindMVA("sun")); TMVA::MethodFisher * KLCoeffsThermal = dynamic_cast(readerKLCoeffs.FindMVA("thermal")); // Fill response tree. for (int entryNum = 0; entryNum < numEntries; ++entryNum) { sampleSumTree -> GetEntry(entryNum); sampleKLCoeffsTree -> GetEntry(entryNum); eventNumberInt = int(eventNumber); cohKLMagIdx = int(coh.KLMagIdx); deconvKLMagIdx = int(deconv.KLMagIdx); bothKLMagIdx = int(both.KLMagIdx); deconvImpulsivity = float(sampleSum -> mostImpulsiveDeconvolvedFiltered(2).impulsivityMeasure); cohImpulsivity = float(sampleSum -> mostImpulsiveCoherentFiltered(2).impulsivityMeasure); deconvLinPolFrac = float(sampleSum -> mostImpulsiveDeconvolvedFiltered(2).linearPolFrac()); cohLinPolFrac = float(sampleSum -> mostImpulsiveCoherentFiltered(2).linearPolFrac()); impulsivityDiff = float(pow(sampleSum -> mostImpulsiveDeconvolvedFiltered(2).impulsivityMeasure, 2) - pow(sampleSum -> mostImpulsiveCoherentFiltered(2).impulsivityMeasure, 2)); peakHilbertDiff = float(sampleSum -> mostImpulsiveDeconvolvedFiltered(2).peakHilbert - sampleSum -> mostImpulsiveCoherentFiltered(2).peakHilbert); deconvInstLinPolFrac = float(sampleSum -> mostImpulsiveDeconvolvedFiltered(2).instantaneousLinearPolFrac()); cohInstLinPolFrac = float(sampleSum -> mostImpulsiveCoherentFiltered(2).instantaneousLinearPolFrac()); deconvPWG = float(sampleSum -> mostImpulsiveDeconvolvedFiltered(2).fracPowerWindowGradient()); cohKLMag = float(coh.KLMag[coh.KLMagIdx]); deconvKLMag = float(deconv.KLMag[deconv.KLMagIdx]); bothKLMag = float(both.KLMag[both.KLMagIdx]); // Request responses. Peng.WAISHPol = readerPeng.EvaluateMVA("WAISHPol"); Peng.WAISVPol = readerPeng.EvaluateMVA("WAISVPol"); Peng.HiCal2A = readerPeng.EvaluateMVA("HiCal2A"); Peng.HiCal2B = readerPeng.EvaluateMVA("HiCal2B"); Peng.iceMC = readerPeng.EvaluateMVA("iceMC"); Peng.RCP = readerPeng.EvaluateMVA("RCP"); Peng.LCP = readerPeng.EvaluateMVA("LCP"); Peng.payloadBlast = readerPeng.EvaluateMVA("payloadBlast"); Peng.sun = readerPeng.EvaluateMVA("sun"); Peng.thermal = readerPeng.EvaluateMVA("thermal"); Andrew.WAISHPol = readerAndrew.EvaluateMVA("WAISHPol"); Andrew.WAISVPol = readerAndrew.EvaluateMVA("WAISVPol"); Andrew.HiCal2A = readerAndrew.EvaluateMVA("HiCal2A"); Andrew.HiCal2B = readerAndrew.EvaluateMVA("HiCal2B"); Andrew.iceMC = readerAndrew.EvaluateMVA("iceMC"); Andrew.RCP = readerAndrew.EvaluateMVA("RCP"); Andrew.LCP = readerAndrew.EvaluateMVA("LCP"); Andrew.payloadBlast = readerAndrew.EvaluateMVA("payloadBlast"); Andrew.sun = readerAndrew.EvaluateMVA("sun"); Andrew.thermal = readerAndrew.EvaluateMVA("thermal"); KLCoeffs.WAISHPol = readerKLCoeffs.EvaluateMVA("WAISHPol"); KLCoeffs.WAISVPol = readerKLCoeffs.EvaluateMVA("WAISVPol"); KLCoeffs.HiCal2A = readerKLCoeffs.EvaluateMVA("HiCal2A"); KLCoeffs.HiCal2B = readerKLCoeffs.EvaluateMVA("HiCal2B"); KLCoeffs.iceMC = readerKLCoeffs.EvaluateMVA("iceMC"); KLCoeffs.RCP = readerKLCoeffs.EvaluateMVA("RCP"); KLCoeffs.LCP = readerKLCoeffs.EvaluateMVA("LCP"); KLCoeffs.payloadBlast = readerKLCoeffs.EvaluateMVA("payloadBlast"); KLCoeffs.sun = readerKLCoeffs.EvaluateMVA("sun"); KLCoeffs.thermal = readerKLCoeffs.EvaluateMVA("thermal"); PengCDF.WAISHPol = PengWAISHPol -> GetRarity(Peng.WAISHPol, TMVA::Types::kSignal); PengCDF.WAISVPol = PengWAISVPol -> GetRarity(Peng.WAISVPol, TMVA::Types::kSignal); PengCDF.HiCal2A = PengHiCal2A -> GetRarity(Peng.HiCal2A, TMVA::Types::kSignal); PengCDF.HiCal2B = PengHiCal2B -> GetRarity(Peng.HiCal2B, TMVA::Types::kSignal); PengCDF.iceMC = PengIceMC -> GetRarity(Peng.iceMC, TMVA::Types::kSignal); PengCDF.RCP = PengRCP -> GetRarity(Peng.RCP, TMVA::Types::kSignal); PengCDF.LCP = PengLCP -> GetRarity(Peng.LCP, TMVA::Types::kSignal); PengCDF.payloadBlast = PengPayloadBlast -> GetRarity(Peng.payloadBlast, TMVA::Types::kSignal); PengCDF.sun = PengSun -> GetRarity(Peng.sun, TMVA::Types::kSignal); PengCDF.thermal = PengThermal -> GetRarity(Peng.thermal, TMVA::Types::kSignal); AndrewCDF.WAISHPol = AndrewWAISHPol -> GetRarity(Andrew.WAISHPol, TMVA::Types::kSignal); AndrewCDF.WAISVPol = AndrewWAISVPol -> GetRarity(Andrew.WAISVPol, TMVA::Types::kSignal); AndrewCDF.HiCal2A = AndrewHiCal2A -> GetRarity(Andrew.HiCal2A, TMVA::Types::kSignal); AndrewCDF.HiCal2B = AndrewHiCal2B -> GetRarity(Andrew.HiCal2B, TMVA::Types::kSignal); AndrewCDF.iceMC = AndrewIceMC -> GetRarity(Andrew.iceMC, TMVA::Types::kSignal); AndrewCDF.RCP = AndrewRCP -> GetRarity(Andrew.RCP, TMVA::Types::kSignal); AndrewCDF.LCP = AndrewLCP -> GetRarity(Andrew.LCP, TMVA::Types::kSignal); AndrewCDF.payloadBlast = AndrewPayloadBlast -> GetRarity(Andrew.payloadBlast, TMVA::Types::kSignal); AndrewCDF.sun = AndrewSun -> GetRarity(Andrew.sun, TMVA::Types::kSignal); AndrewCDF.thermal = AndrewThermal -> GetRarity(Andrew.thermal, TMVA::Types::kSignal); KLCoeffsCDF.WAISHPol = KLCoeffsWAISHPol -> GetRarity(KLCoeffs.WAISHPol, TMVA::Types::kSignal); KLCoeffsCDF.WAISVPol = KLCoeffsWAISVPol -> GetRarity(KLCoeffs.WAISVPol, TMVA::Types::kSignal); KLCoeffsCDF.HiCal2A = KLCoeffsHiCal2A -> GetRarity(KLCoeffs.HiCal2A, TMVA::Types::kSignal); KLCoeffsCDF.HiCal2B = KLCoeffsHiCal2B -> GetRarity(KLCoeffs.HiCal2B, TMVA::Types::kSignal); KLCoeffsCDF.iceMC = KLCoeffsIceMC -> GetRarity(KLCoeffs.iceMC, TMVA::Types::kSignal); KLCoeffsCDF.RCP = KLCoeffsRCP -> GetRarity(KLCoeffs.RCP, TMVA::Types::kSignal); KLCoeffsCDF.LCP = KLCoeffsLCP -> GetRarity(KLCoeffs.LCP, TMVA::Types::kSignal); KLCoeffsCDF.payloadBlast = KLCoeffsPayloadBlast -> GetRarity(KLCoeffs.payloadBlast, TMVA::Types::kSignal); KLCoeffsCDF.sun = KLCoeffsSun -> GetRarity(KLCoeffs.sun, TMVA::Types::kSignal); KLCoeffsCDF.thermal = KLCoeffsThermal -> GetRarity(KLCoeffs.thermal, TMVA::Types::kSignal); responseTree.Fill(); } // Write responseTree to file, then close it. sampleResponseFile.cd(); // responseTree.BuildIndex("run", "eventNumber"); responseTree.Write(); sampleResponseFile.Close(); sampleKLCoeffsFile.Close(); sampleSumFile.Close(); }