-
Notifications
You must be signed in to change notification settings - Fork 1
/
michelEventLoop.txt
1438 lines (1267 loc) · 86.1 KB
/
michelEventLoop.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#define MC_OUT_FILE_NAME "runEventLoopMC_July202022_EavailComp_q3bin6_reweight.root"
#define DATA_OUT_FILE_NAME "runEventLoopData_July202022_EavailComp_q3bin6_reweight.root"
#define MC_SIDE_FILE_NAME "runEventLoopMC_July202022_Sideband_q3bin6_reweight.root"
#define DATA_SIDE_FILE_NAME "runEventLoopDATA_July202022_Sideband_q3bin6_reweight.root"
#define USAGE \
"\n*** USAGE ***\n"\
"runEventLoop <dataPlaylist.txt> <mcPlaylist.txt>\n\n"\
"*** Explanation ***\n"\
"Reduce MasterAnaDev AnaTuples to event selection histograms to extract a\n"\
"single-differential inclusive cross section for the 2021 MINERvA 101 tutorial.\n\n"\
"*** The Input Files ***\n"\
"Playlist files are plaintext files with 1 file name per line. Filenames may be\n"\
"xrootd URLs or refer to the local filesystem. The first playlist file's\n"\
"entries will be treated like data, and the second playlist's entries must\n"\
"have the \"Truth\" tree to use for calculating the efficiency denominator.\n\n"\
"*** Output ***\n"\
"Produces a two files, " MC_OUT_FILE_NAME " and " DATA_OUT_FILE_NAME ", with\n"\
"all histograms needed for the ExtractCrossSection program also built by this\n"\
"package. You'll need a .rootlogon.C that loads ROOT object definitions from\n"\
"PlotUtils to access systematics information from these files.\n\n"\
"*** Environment Variables ***\n"\
"Setting up this package appends to PATH and LD_LIBRARY_PATH. PLOTUTILSROOT,\n"\
"MPARAMFILESROOT, and MPARAMFILES must be set according to the setup scripts in\n"\
"those packages for systematics and flux reweighters to function.\n"\
"If MNV101_SKIP_SYST is defined at all, output histograms will have no error bands.\n"\
"This is useful for debugging the CV and running warping studies.\n\n"\
"*** Return Codes ***\n"\
"0 indicates success. All histograms are valid only in this case. Any other\n"\
"return code indicates that histograms should not be used. Error messages\n"\
"about what went wrong will be printed to stderr. So, they'll end up in your\n"\
"terminal, but you can separate them from everything else with something like:\n"\
"\"runEventLoop data.txt mc.txt 2> errors.txt\"\n"
enum ErrorCodes
{
success = 0,
badCmdLine = 1,
badInputFile = 2,
badFileRead = 3,
badOutputFile = 4
};
//PlotUtils includes
//No junk from PlotUtils please! I already
//know that MnvH1D does horrible horrible things.
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Woverloaded-virtual"
//Includes from this package
#include "event/CVUniverse.h"
#include "event/MichelEvent.h"
#include "systematics/Systematics.h"
#include "cuts/MaxPzMu.h"
#include "util/Variable.h"
#include "util/Variable2D.h"
#include "util/GetFluxIntegral.h"
#include "util/GetPlaylist.h"
#include "cuts/SignalDefinition.h"
//#include "util/TruthInteractionStudies.h"
#include "util/PionFSCategory.h"
#include "cuts/q3RecoCut.h"
#include "studies/Study.h"
//#include "studies/PerMichelVarByGENIELabel.h"
#include "studies/PerMichelEventVarByGENIELabel.h"
#include "studies/PerMichel2DVar.h"
#include "studies/PerMichel2DVarbin.h"
#include "studies/PerMichelVarVec.h"
#include "studies/PerMichelEvent2DVarBin.h"
#include "studies/PerEventVarBin.h"
#include "studies/PerMichelVarVecFSPart.h"
#include "cuts/hasMichel.h"
#include "event/Michel.h"
#include "cuts/BestMichelDistance2D.h"
#include "cuts/VtxMatchFirst.h"
#include "cuts/hasTruePion.h"
#include "cuts/PTRangeReco.h"
#include "cuts/GetClosestMichel.h"
#include "cuts/Distance2DSideband.h"
#include "event/SetDistanceMichelSideband.h"
#include "event/SetDistanceMichelSelection.h"
#include "event/GetClosestMichel.h"
//#include "Binning.h" //TODO: Fix me
//PlotUtils includes
#include "PlotUtils/makeChainWrapper.h"
#include "PlotUtils/HistWrapper.h"
#include "PlotUtils/Hist2DWrapper.h"
#include "PlotUtils/MacroUtil.h"
#include "PlotUtils/MnvPlotter.h"
#include "PlotUtils/CCInclusiveCuts.h"
#include "PlotUtils/CCInclusiveSignal.h"
#include "PlotUtils/CrashOnROOTMessage.h" //Sets up ROOT's debug callbacks by itself
#include "PlotUtils/Cutter.h"
#include "PlotUtils/Model.h"
#include "PlotUtils/FluxAndCVReweighter.h"
#include "PlotUtils/GENIEReweighter.h"
#include "PlotUtils/LowRecoil2p2hReweighter.h"
#include "PlotUtils/RPAReweighter.h"
#include "PlotUtils/MINOSEfficiencyReweighter.h"
#include "PlotUtils/TargetUtils.h"
#include "../MAT-MINERvA/weighters/PionReweighter.h"
#include "PlotUtils/LowRecPionSignal.h"
#pragma GCC diagnostic pop
//ROOT includes
#include "TParameter.h"
//c++ includes
#include <iostream>
#include <cstdlib> //getenv()
//==============================================================================
// Loop and Fill
//==============================================================================
void LoopAndFillEventSelection(
PlotUtils::ChainWrapper* chain,
std::map<std::string, std::vector<CVUniverse*> > error_bands,
std::vector<Variable*> vars,
std::vector<Variable2D*> vars2D,
std::vector<Study*> studies,
std::vector<Variable*> sidevars,
std::vector<Variable2D*> sidevars2D,
std::vector<Study*> sideband_studies,
PlotUtils::Cutter<CVUniverse, MichelEvent>& michelcuts,
PlotUtils::Model<CVUniverse, MichelEvent>& model)
{
assert(!error_bands["cv"].empty() && "\"cv\" error band is empty! Can't set Model weight.");
auto& cvUniv = error_bands["cv"].front();
std::cout << "Starting MC reco loop...\n";
const int nEntries = chain->GetEntries(); // TODO: July 10 CHANGE BACK TO GETENTRIES
for (int i=0; i < nEntries; ++i)
{
if(i%1000==0) std::cout << i << " / " << nEntries << "\r" << std::flush;
//std::cout << "Now Printing for Event " << i << std::endl;
MichelEvent cvEvent;
cvUniv->SetEntry(i);
model.SetEntry(*cvUniv, cvEvent);
const double cvWeight =model.GetWeight(*cvUniv, cvEvent);// TODO: Put this model weight back. model.GetWeight(*cvUniv, cvEvent);
//=========================================
// Systematics loop(s)
//=========================================
for (auto band : error_bands)
{
std::vector<CVUniverse*> error_band_universes = band.second;
for (auto universe : error_band_universes)
{
//if (universe->ShortName() != "cv") continue;
MichelEvent myevent; // make sure your event is inside the error band loop.
//if (universe->ShortName() != "cv") continue;
// Tell the Event which entry in the TChain it's looking at
universe->SetEntry(i);
const double weight2 = model.GetWeight(*universe, myevent);
const auto cutResults = michelcuts.isMCSelected(*universe, myevent, cvWeight);
//const auto cutResults = michelcuts.isDataSelected(*universe, myevent);
//if (universe->ShortName() != "cv") continue;
if (!cutResults.all()) continue;
if (cutResults.all()){
setDistanceMichelSelection(*universe, myevent, 150.);
setClosestMichel(*universe, myevent,0);
if (!myevent.m_nmichelspass.empty()){
for(auto& var: vars) {
//std::cout << "Filling Var for this universe " << i << std::endl;
//std::cout << " Var Value is " << var->GetRecoValue(*universe) << std::endl;
var->selectedMCReco->FillUniverse(universe, var->GetRecoValue(*universe), weight2); //"Fake data" for closure
(*var->m_MChists)[universe->GetInteractionType()].FillUniverse(universe, var->GetRecoValue(*universe), weight2);
var->FillCategHistos(*universe,var->GetRecoValue(*universe), weight2);
//var->migration->FillUniverse(universe, var->GetRecoValue(*universe), var->GetTrueValue(*universe), weight2);
var->efficiencyNumerator->FillUniverse(universe, var->GetTrueValue(*universe), weight2);
}
for(auto& var: vars2D) {
var->mcTotalHist->FillUniverse(universe, var->GetRecoValueX(*universe), var->GetRecoValueY(*universe), weight2);
(*var->m_MChists)[universe->GetInteractionType()].FillUniverse(universe, var->GetRecoValueX(*universe), var->GetRecoValueY(*universe), weight2);
var->FillCategHistos(*universe,var->GetRecoValueX(*universe), var->GetRecoValueY(*universe), weight2);
//var->efficiencyNumerator->FillUniverse(universe, var->GetTrueValueX(*universe), var->GetTrueValueY(*universe), weight2);
}
for(auto& study: studies) study->SelectedSignal(*universe, myevent, weight2);
const bool isSignal = michelcuts.isSignal(*universe, weight2);
if(isSignal){
for(auto& var: vars)
{
var->migration->FillUniverse(universe, var->GetRecoValue(*universe), var->GetTrueValue(*universe), weight2);
var->efficiencyNumerator->FillUniverse(universe, var->GetTrueValue(*universe), weight2);
var->selectedSignalReco->FillUniverse(universe, var->GetRecoValue(*universe), weight2); //Efficiency numerator in reco variables. Useful for warping studies
}
for(auto& var: vars2D)
{
var->efficiencyNumerator->FillUniverse(universe, var->GetTrueValueX(*universe), var->GetTrueValueY(*universe), weight2);
}
} // end of if Signal()
else{
int bkgd_ID = -1;
if (universe->GetCurrent()==2)bkgd_ID=0;
else bkgd_ID=1;
for(auto& var: vars){
(*var->m_backgroundHists)[bkgd_ID].FillUniverse(universe, var->GetRecoValue(*universe), weight2);
}
for(auto& var: vars2D) (*var->m_backgroundHists)[bkgd_ID].FillUniverse(universe, var->GetRecoValueX(*universe), var->GetRecoValueY(*universe), weight2);
} // End of else statement for if Signal
} // End of if SelectedMichels Not Empty
else{
//myevent.m_nmichels = myevent.m_sidebandpass;
setDistanceMichelSidebands(*universe, myevent, 150., 500.);
setClosestMichel(*universe, myevent,1);
if (!myevent.m_sidebandpass.empty()){
//const double weight2 = model.GetWeight(*universe, myevent);
for(auto& var: sidevars) {
var->selectedMCReco->FillUniverse(universe, var->GetRecoValue(*universe), weight2); //"Fake data" for closure
(*var->m_MChists)[universe->GetInteractionType()].FillUniverse(universe, var->GetRecoValue(*universe), weight2);
var->FillCategHistos(*universe,var->GetRecoValue(*universe), weight2);
var->migration->FillUniverse(universe, var->GetRecoValue(*universe), var->GetTrueValue(*universe), weight2);
var->efficiencyNumerator->FillUniverse(universe, var->GetTrueValue(*universe), weight2);
}
for(auto& var: sidevars2D) {
var->mcTotalHist->FillUniverse(universe, var->GetRecoValueX(*universe), var->GetRecoValueY(*universe), weight2);
(*var->m_MChists)[universe->GetInteractionType()].FillUniverse(universe, var->GetRecoValueX(*universe), var->GetRecoValueY(*universe), weight2);
var->FillCategHistos(*universe,var->GetRecoValueX(*universe), var->GetRecoValueY(*universe), weight2);
var->efficiencyNumerator->FillUniverse(universe, var->GetTrueValueX(*universe), var->GetTrueValueY(*universe), weight2);
}
for(auto& study: sideband_studies) study->SelectedSignal(*universe, myevent, weight2);
}
} // end of else if (!cutResults[0] && evt.sideband == 1) // To fill Sideband Variables
} // If event passes PreCuts
} // End band's universe loop
} // End Band loop
} //End entries loop
std::cout << "Finished MC reco loop.\n";
}
void LoopAndFillData( PlotUtils::ChainWrapper* data,
std::vector<CVUniverse*> data_band,
std::vector<Variable*> vars,
std::vector<Variable2D*> vars2D,
std::vector<Study*> studies,
std::vector<Variable*> sidevars,
std::vector<Variable2D*> sidevars2D,
std::vector<Study*> data_sidebands,
PlotUtils::Cutter<CVUniverse, MichelEvent>& michelcuts)
{
std::cout << "Starting data loop...\n";
const int nEntries = data->GetEntries(); // TODO: July 10 CHANGE BACK TO GEtENTRIES
for (int i=0; i <nEntries; ++i) {
//std::cout << "Now Printing for Event " << i << std::endl;
//for (auto universe : data_band) {
const auto universe = data_band.front();
universe->SetEntry(i);
if(i%1000==0) std::cout << i << " / " << nEntries << "\r" << std::flush;
//std::cout << "Creating Michel Event" << std::endl;
//if (universe->ShortName() != "cv") continue;
MichelEvent myevent;
//std::cout << "Applying Cuts to Data Event " << std::endl;
const auto cutResults = michelcuts.isDataSelected(*universe, myevent);
if (!cutResults.all()) continue;
if (cutResults.all()){
setDistanceMichelSelection(*universe, myevent, 150.);
setClosestMichel(*universe, myevent,0);
if (!myevent.m_nmichelspass.empty()){
//std::cout << "Filling Data STudies" << std::endl;
for(auto& study: studies) study->Selected(*universe, myevent, 1);
for(auto& var: vars)
{
//std::cout << "Filling Data Var for This universe " << i << std::endl;
//std::cout << " Var Value is " << var->GetRecoValue(*universe) << std::endl;
var->dataHist->FillUniverse(universe, var->GetRecoValue(*universe), 1);
}
for(auto& var: vars2D)
{
var->dataHist->FillUniverse(universe, var->GetRecoValueX(*universe), var->GetRecoValueY(*universe), 1);
}
} // End of if (michelcuts.isDataSelected(*universe, myevent).all()) -> filling Selection
else{
setDistanceMichelSidebands(*universe, myevent, 150., 500.);
setClosestMichel(*universe, myevent,1);
if (!myevent.m_sidebandpass.empty()){
for(auto& study: data_sidebands) study->Selected(*universe, myevent, 1);
for(auto& var: sidevars)
{
var->dataHist->FillUniverse(universe, var->GetRecoValue(*universe), 1);
}
for(auto& var: sidevars2D)
{
var->dataHist->FillUniverse(universe, var->GetRecoValueX(*universe), var->GetRecoValueY(*universe), 1);
}
} // End of else for Filling Sideband Variables
} // End of else if no michelspass
}// End of PreCuts Pass
//} // End of Data Band
} // End of Entries
std::cout << "Finished data loop.\n";
}
void LoopAndFillEffDenom( PlotUtils::ChainWrapper* truth,
std::map<std::string, std::vector<CVUniverse*> > truth_bands,
std::vector<Variable*> vars,
std::vector<Variable2D*> vars2D,
std::vector<Study*> studies,
PlotUtils::Cutter<CVUniverse, MichelEvent>& michelcuts,
PlotUtils::Model<CVUniverse, MichelEvent>& model)
{
assert(!truth_bands["cv"].empty() && "\"cv\" error band is empty! Could not set Model entry.");
auto& cvUniv = truth_bands["cv"].front();
std::cout << "Starting efficiency denominator loop...\n";
const int nEntries = truth->GetEntries();
for (int i=0; i<nEntries; ++i)
{
if(i%1000==0) std::cout << i << " / " << nEntries << "\r" << std::flush;
MichelEvent cvEvent;
cvUniv->SetEntry(i);
model.SetEntry(*cvUniv, cvEvent);
const double cvWeight = model.GetWeight(*cvUniv, cvEvent);
//=========================================
// Systematics loop(s)
//=========================================
for (auto band : truth_bands)
{
std::vector<CVUniverse*> truth_band_universes = band.second;
for (auto universe : truth_band_universes)
{
MichelEvent myevent; //Only used to keep the Model happy
// Tell the Event which entry in the TChain it's looking at
universe->SetEntry(i);
if (!michelcuts.isEfficiencyDenom(*universe, cvWeight)) continue; //Weight is ignored for isEfficiencyDenom() in all but the CV universe
const double weight = model.GetWeight(*universe, myevent); //Only calculate the weight for events that will use it
//Fill efficiency denominator now:
for(auto var: vars)
{
var->efficiencyDenominator->FillUniverse(universe, var->GetTrueValue(*universe), weight);
}
// Fill Studies denominator:
// for(auto& study: studies) study->SelectedSignal(*universe, cvEvent, weight);
for(auto var: vars2D)
{
var->efficiencyDenominator->FillUniverse(universe, var->GetTrueValueX(*universe), var->GetTrueValueY(*universe), weight);
}
}
}
}
std::cout << "Finished efficiency denominator loop.\n";
}
//Returns false if recoTreeName could not be inferred
bool inferRecoTreeNameAndCheckTreeNames(const std::string& mcPlaylistName, const std::string& dataPlaylistName, std::string& recoTreeName)
{
const std::vector<std::string> knownTreeNames = {"Truth", "Meta", "Header"};
bool areFilesOK = false;
std::ifstream playlist(mcPlaylistName);
std::string firstFile = "";
playlist >> firstFile;
auto testFile = TFile::Open(firstFile.c_str());
if(!testFile)
{
std::cerr << "Failed to open the first MC file at " << firstFile << "\n";
return false;
}
//Does the MC playlist have the Truth tree? This is needed for the efficiency denominator.
const auto truthTree = testFile->Get("Truth");
if(truthTree == nullptr || !truthTree->IsA()->InheritsFrom(TClass::GetClass("TTree")))
{
std::cerr << "Could not find the \"Truth\" tree in MC file named " << firstFile << "\n";
return false;
}
//Figure out what the reco tree name is
for(auto key: *testFile->GetListOfKeys())
{
if(static_cast<TKey*>(key)->ReadObj()->IsA()->InheritsFrom(TClass::GetClass("TTree"))
&& std::find(knownTreeNames.begin(), knownTreeNames.end(), key->GetName()) == knownTreeNames.end())
{
recoTreeName = key->GetName();
areFilesOK = true;
}
}
delete testFile;
testFile = nullptr;
//Make sure the data playlist's first file has the same reco tree
playlist.open(dataPlaylistName);
playlist >> firstFile;
testFile = TFile::Open(firstFile.c_str());
if(!testFile)
{
std::cerr << "Failed to open the first data file at " << firstFile << "\n";
return false;
}
const auto recoTree = testFile->Get(recoTreeName.c_str());
if(recoTree == nullptr || !recoTree->IsA()->InheritsFrom(TClass::GetClass("TTree")))
{
std::cerr << "Could not find the \"" << recoTreeName << "\" tree in data file named " << firstFile << "\n";
return false;
}
return areFilesOK;
}
//==============================================================================
// Main
//==============================================================================
int main(const int argc, const char** argv)
{
TH1::AddDirectory(false);
//Validate input.
//I expect a data playlist file name and an MC playlist file name which is exactly 2 arguments.
const int nArgsExpected = 2;
if(argc != nArgsExpected + 1) //argc is the size of argv. I check for number of arguments + 1 because
//argv[0] is always the path to the executable.
{
std::cerr << "Expected " << nArgsExpected << " arguments, but got " << argc - 1 << "\n" << USAGE << "\n";
return badCmdLine;
}
//One playlist must contain only MC files, and the other must contain only data files.
//Only checking the first file in each playlist because opening each file an extra time
//remotely (e.g. through xrootd) can get expensive.
//TODO: Look in INSTALL_DIR if files not found?
const std::string mc_file_list = argv[2],
data_file_list = argv[1];
//Check that necessary TTrees exist in the first file of mc_file_list and data_file_list
std::string reco_tree_name;
if(!inferRecoTreeNameAndCheckTreeNames(mc_file_list, data_file_list, reco_tree_name))
{
std::cerr << "Failed to find required trees in MC playlist " << mc_file_list << " and/or data playlist " << data_file_list << ".\n" << USAGE << "\n";
return badInputFile;
}
const bool doCCQENuValidation = (reco_tree_name == "CCQENu"); //Enables extra histograms and might influence which systematics I use.
const bool is_grid = false;
PlotUtils::MacroUtil options(reco_tree_name, mc_file_list, data_file_list, "minervame1A", true);
std::cout << options.m_mc->GetChain()->GetName() << std::endl;
options.m_plist_string = util::GetPlaylist(*options.m_mc, true); //TODO: Put GetPlaylist into PlotUtils::MacroUtil
// You're required to make some decisions
PlotUtils::MinervaUniverse::SetNuEConstraint(true);
PlotUtils::MinervaUniverse::SetPlaylist(options.m_plist_string); //TODO: Infer this from the files somehow?
PlotUtils::MinervaUniverse::SetAnalysisNuPDG(14);
PlotUtils::MinervaUniverse::SetNFluxUniverses(100);
PlotUtils::MinervaUniverse::SetZExpansionFaReweight(false);
//Now that we've defined what a cross section is, decide which sample and model
//we're extracting a cross section for.
PlotUtils::Cutter<CVUniverse, MichelEvent>::reco_t preCuts;
PlotUtils::Cutter<CVUniverse, MichelEvent>::reco_t sidebands;
PlotUtils::Cutter<CVUniverse, MichelEvent>::truth_t signalDefinition, phaseSpace;
//const double minZ = 5980, maxZ = 8590.07, apothem = 850.; // All in mm
const double minZ = 5980, maxZ = 8422, apothem = 850; //All in mm
//preCuts.emplace_back(new reco::ZRange<CVUniverse, MichelEvent>("Tracker", minZ, maxZ));
preCuts.emplace_back(new reco::Apothem<CVUniverse, MichelEvent>(apothem));
preCuts.emplace_back(new reco::ZRange<CVUniverse, MichelEvent>("Tracker", minZ, maxZ));
preCuts.emplace_back(new reco::MaxMuonAngle<CVUniverse, MichelEvent>(20.));
preCuts.emplace_back(new reco::HasMINOSMatch<CVUniverse, MichelEvent>());
preCuts.emplace_back(new reco::NoDeadtime<CVUniverse, MichelEvent>(1, "Deadtime"));
preCuts.emplace_back(new reco::IsNeutrino<CVUniverse, MichelEvent>());
//preCuts.emplace_back(new Q3RangeReco<CVUniverse, MichelEvent>(0.0,1.2));
preCuts.emplace_back(new PTRangeReco<CVUniverse, MichelEvent>(0.0,1.0));
//preCuts.emplace_back(new hasMichel<CVUniverse, MichelEvent>());
//preCuts.emplace_back(new BestMichelDistance2D<CVUniverse, MichelEvent>(100.));
//preCuts.emplace_back(new VtxMatchFirst<CVUniverse, MichelEvent>(200., 102.));
//preCuts.emplace_back(new NPiCut<CVUniverse, MichelEvent>(1));
preCuts.emplace_back(new hasMichel<CVUniverse, MichelEvent>());
//preCuts.emplace_back(new Distance2DSideband<CVUniverse, MichelEvent>(500.));
//preCuts.emplace_back(new BestMichelDistance2D<CVUniverse, MichelEvent>(150.));
//preCuts.emplace_back(new GetClosestMichel<CVUniverse, MichelEvent>(0));
//nosidebands.emplace_back(new BestMichelDistance2D<CVUniverse, MichelEvent>(150.));
//nosidebands.emplace_back(new GetClosestMichel<CVUniverse, MichelEvent>(0));
TFile* mc_MichelStudies = TFile::Open("July202022_EavailComp_noreweight_q3bin6_MC.root", "RECREATE");
TFile* data_MichelStudies = TFile::Open("July202022_EavailComp_noreweight_q3bin6_data.root", "RECREATE");
TFile* mc_SidebandStudies = TFile::Open("July202022_Sideband_noreweight_q3bin6_MC.root", "RECREATE");
TFile* data_SidebandStudies = TFile::Open("July202022_Sideband_noreweight_q3bin6_data.root", "RECREATE");
signalDefinition.emplace_back(new truth::IsNeutrino<CVUniverse>());
signalDefinition.emplace_back(new truth::IsCC<CVUniverse>());
signalDefinition.emplace_back(new truth::HasPion<CVUniverse>());
//signalDefinition.emplace_back(new Q3Limit<CVUniverse>(0.0, 1.20));
//signalDefinition.emplace_back(new hasTruePion<CVUniverse>());
phaseSpace.emplace_back(new truth::ZRange<CVUniverse>("Tracker", minZ, maxZ));
phaseSpace.emplace_back(new truth::Apothem<CVUniverse>(apothem));
phaseSpace.emplace_back(new truth::MuonAngle<CVUniverse>(20.));
phaseSpace.emplace_back(new truth::PZMuMin<CVUniverse>(1500.));
phaseSpace.emplace_back(new truth::pTRangeLimit<CVUniverse>(0., 1.0));
PlotUtils::Cutter<CVUniverse, MichelEvent> mycuts(std::move(preCuts), std::move(sidebands) , std::move(signalDefinition),std::move(phaseSpace));
std::vector<std::unique_ptr<PlotUtils::Reweighter<CVUniverse, MichelEvent>>> MnvTunev1;
MnvTunev1.emplace_back(new PlotUtils::FluxAndCVReweighter<CVUniverse, MichelEvent>());
MnvTunev1.emplace_back(new PlotUtils::GENIEReweighter<CVUniverse, MichelEvent>(true, false));
MnvTunev1.emplace_back(new PlotUtils::LowRecoil2p2hReweighter<CVUniverse, MichelEvent>());
MnvTunev1.emplace_back(new PlotUtils::MINOSEfficiencyReweighter<CVUniverse, MichelEvent>());
MnvTunev1.emplace_back(new PlotUtils::RPAReweighter<CVUniverse, MichelEvent>());
//TODO: Add my pion reweighter here. - Mehreen S. Nov 22, 2021
//MnvTunev1.emplace_back(new PlotUtils::PionReweighter<CVUniverse,MichelEvent>());
PlotUtils::Model<CVUniverse, MichelEvent> model(std::move(MnvTunev1));
// Make a map of systematic universes
// Leave out systematics when making validation histograms
const bool doSystematics = (getenv("MNV101_SKIP_SYST") == nullptr);
if(!doSystematics){
std::cout << "Skipping systematics (except 1 flux universe) because environment variable MNV101_SKIP_SYST is set.\n";
PlotUtils::MinervaUniverse::SetNFluxUniverses(2); //Necessary to get Flux integral later... Doesn't work with just 1 flux universe though because _that_ triggers "spread errors".
}
std::map< std::string, std::vector<CVUniverse*> > error_bands;
if(doSystematics) error_bands = GetStandardSystematics(options.m_mc);
else{
std::map<std::string, std::vector<CVUniverse*> > band_flux = PlotUtils::GetFluxSystematicsMap<CVUniverse>(options.m_mc, CVUniverse::GetNFluxUniverses());
error_bands.insert(band_flux.begin(), band_flux.end()); //Necessary to get flux integral later...
}
error_bands["cv"] = {new CVUniverse(options.m_mc)};
std::map< std::string, std::vector<CVUniverse*> > truth_bands;
if(doSystematics) truth_bands = GetStandardSystematics(options.m_truth);
truth_bands["cv"] = {new CVUniverse(options.m_truth)};
std::vector<double> dansPTBins = {0, 0.075, 0.10, 0.15, 0.20, 0.30, 0.4, 0.50,0.60 , 0.7, 0.80,0.9, 1.,1.1, 1.2, 1.3, 1.4, 1.5},
dansPzBins = {1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 6, 7, 8, 9, 10, 15, 20, 40, 60},
robsEmuBins = {0,1,2,3,4,5,7,9,12,15,18,22,36,50,75,80},
mehreenQ3Bins = {0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4},
robsRecoilBins;
int nq3mbins = mehreenQ3Bins.size() -1;
std::vector<double> tpibins = {0, 4., 8., 12., 16., 20., 24., 28., 32., 36., 40., 46., 52.,60., 70., 80., 100., 125.,150., 175., 200., 225., 250., 275., 300., 325., 350., 400., 500., 1000.};
std::vector<double> rangebins = {0, 4., 8., 12., 16., 20., 24., 28., 32., 36., 40., 44., 50., 56., 62., 70., 80.,90., 100., 110., 120., 140., 160., 180., 200., 220., 240., 260., 280., 300., 325., 350., 375., 400., 450., 500., 550., 600., 650., 700., 800., 900., 1000., 1200., 1400., 1800., 2400.};
//std::vector<double> tpibins = {0., 4., 8., 12., 16., 20., 24., 28., 32., 36., 40., 46., 52., 60., 70., 80., 90., 100., 120., 140., 160., 180., 200., 220., 240., 260.,280., 300., 320., 340., 360., 380., 400., 500., 1000.};
//std::vector<double> rangebins = {0., 4., 8., 12., 16., 20., 24., 28., 32., 36., 40., 46., 52., 60., 70., 80., 90., 100., 120., 140., 160., 180., 200., 220., 260., 280., 300., 320., 340., 360., 380., 400., 425., 450., 475., 500., 600., 700., 800., 1000., 1200., 1400., 1600., 2000., 2400.};
std::vector<double> recoilbins = {0.0, 150., 200., 300., 400., 500., 600., 800., 1000., 1200., 1400., 1600.};
const double robsRecoilBinWidth = 50; //MeV
for(int whichBin = 0; whichBin < 30 + 1; ++whichBin) robsRecoilBins.push_back(robsRecoilBinWidth * whichBin);
std::vector<Variable*> vars = {
new Variable("pTmu", "p_{T, #mu} [GeV/c]", dansPTBins, &CVUniverse::GetMuonPT, &CVUniverse::GetMuonPTTrue),
new Variable("q3", "q3 [GeV]", dansPTBins, &CVUniverse::Getq3, &CVUniverse::GetTrueQ3),
new Variable("q2", "q2 [GeV^2]", dansPTBins, &CVUniverse::GetQ2Reco, &CVUniverse::GetTrueQ2GeV),
new Variable("pzmu", "p_{||, #mu} [GeV/c]", dansPzBins, &CVUniverse::GetMuonPz, &CVUniverse::GetMuonPzTrue),
new Variable("Emu", "E_{#mu} [GeV]", robsEmuBins, &CVUniverse::GetEmuGeV, &CVUniverse::GetElepTrueGeV),
};
std::vector<Variable*> varsagain = {
new Variable("pTmu", "pT [GeV/c]", dansPTBins, &CVUniverse::GetMuonPT, &CVUniverse::GetMuonPTTrue),
new Variable("q3", "q3 [GeV]", mehreenQ3Bins, &CVUniverse::Getq3, &CVUniverse::GetTrueQ3),
new Variable("q2", "q2 [GeV^2]", dansPTBins, &CVUniverse::GetQ2Reco, &CVUniverse::GetTrueQ2GeV),
new Variable("pzmu", "p_z [GeV/c]", dansPzBins, &CVUniverse::GetMuonPz, &CVUniverse::GetMuonPzTrue),
new Variable("Emu", "E_{#mu} [GeV]", robsEmuBins, &CVUniverse::GetEmuGeV, &CVUniverse::GetElepTrueGeV),
new Variable("pTmubins", "pT [GeV/c]", mehreenQ3Bins, &CVUniverse::GetMuonPT, &CVUniverse::GetMuonPTTrue),
};
std::vector<Variable2D*> vars2D;
vars2D.push_back(new Variable2D(*varsagain[6], *varsagain[1])); // ptmu bins vs pt
vars2D.push_back(new Variable2D(*varsagain[6], *varsagain[3]));// ptmubins vs Eavail For some reason this is not printing correctly
//vars2D.push_back(new Variable2D(*varsagain[1], *varsagain[3]));// ptmu vs Eavail
vars2D.push_back(new Variable2D(*varsagain[6], *varsagain[2])); // ptmu bins vs q3
vars2D.push_back(new Variable2D(*varsagain[6], *varsagain[4])); // ptmu bins vs pz
//vars2D.push_back(new Variable2D(*varsagain[1], *varsagain[2])); // ptmu bins vs q2
//vars2D.push_back(new Variable2D(*varsagain[1], *varsagain[4]));
//vars2D.push_back(new Variable2D(*varsagain[1], *varsagain[0]));
//vars2D.push_back(new Variable2D(*varsagain[3], *varsagain[0]));
//vars2D.push_back(new Variable2D(*varsagain[3], *varsagain[3]));
vars2D.push_back(new Variable2D(*varsagain[7], *varsagain[12]));
//vars2D.push_back(new Variable2D(*varsagain[8], *varsagain[13]));
vars2D.push_back(new Variable2D(*varsagain[9], *varsagain[14]));
//vars2D.push_back(new Variable2D(*varsagain[10], *varsagain[15]));
vars2D.push_back(new Variable2D(*varsagain[11], *varsagain[16]));
vars2D.push_back(new Variable2D(*varsagain[3], *varsagain[17])); // Eavail vs true q0 in recoilbins
vars2D.push_back(new Variable2D(*varsagain[18], *varsagain[21])); //New Eavail vs true Eavail in recoil bins y;
vars2D.push_back(new Variable2D(*varsagain[19], *varsagain[22])); // New RecoilE vs true q0 in recoil bins y;
vars2D.push_back(new Variable2D(*varsagain[20], *varsagain[23])); // New Eavail vs true q0 in recoil bins y;
std::vector<Variable*> sidevars = {
new Variable("pTmu", "p_{T, #mu} [GeV/c]", dansPTBins, &CVUniverse::GetMuonPT, &CVUniverse::GetMuonPTTrue),
new Variable("q3", "q3 [GeV]", dansPTBins, &CVUniverse::Getq3, &CVUniverse::GetTrueQ3),
new Variable("q2", "q2 [GeV^2]", dansPTBins, &CVUniverse::GetQ2Reco, &CVUniverse::GetTrueQ2GeV),
new Variable("pzmu", "p_{||, #mu} [GeV/c]", dansPzBins, &CVUniverse::GetMuonPz, &CVUniverse::GetMuonPzTrue),
new Variable("Emu", "E_{#mu} [GeV]", robsEmuBins, &CVUniverse::GetEmuGeV, &CVUniverse::GetElepTrueGeV),
new Variable("q3pTdiff","[GeV]", dansPTBins, &CVUniverse::Recoq3pTdiff, &CVUniverse::GetTrueq3pTdiff),
};
std::vector<Variable*> sidevarsagain = {
new Variable("pTmu", "pT [GeV/c]", dansPTBins, &CVUniverse::GetMuonPT, &CVUniverse::GetMuonPTTrue),
new Variable("q3", "q3 [GeV]", mehreenQ3Bins, &CVUniverse::Getq3, &CVUniverse::GetTrueQ3),
new Variable("q2", "q2 [GeV^2]", dansPTBins, &CVUniverse::GetQ2Reco, &CVUniverse::GetTrueQ2GeV),
new Variable("pzmu", "p_z [GeV/c]", dansPzBins, &CVUniverse::GetMuonPz, &CVUniverse::GetMuonPzTrue),
new Variable("Emu", "E_{#mu} [GeV]", robsEmuBins, &CVUniverse::GetEmuGeV, &CVUniverse::GetElepTrueGeV),
new Variable("pTmubins", "pT [GeV/c]", mehreenQ3Bins, &CVUniverse::GetMuonPT, &CVUniverse::GetMuonPTTrue),
};
std::vector<Variable2D*> sidevars2D;
sidevars2D.push_back(new Variable2D(*sidevarsagain[6], *sidevarsagain[1]));
sidevars2D.push_back(new Variable2D(*sidevarsagain[6], *sidevarsagain[3]));//For some reason this is not printing correctly
sidevars2D.push_back(new Variable2D(*sidevarsagain[1], *sidevarsagain[3]));
sidevars2D.push_back(new Variable2D(*sidevarsagain[6], *sidevarsagain[2]));
sidevars2D.push_back(new Variable2D(*sidevarsagain[6], *sidevarsagain[4]));
sidevars2D.push_back(new Variable2D(*sidevarsagain[1], *sidevarsagain[2]));
sidevars2D.push_back(new Variable2D(*sidevarsagain[1], *sidevarsagain[4]));
sidevars2D.push_back(new Variable2D(*sidevarsagain[1], *sidevarsagain[0]));
sidevars2D.push_back(new Variable2D(*sidevarsagain[3], *sidevarsagain[0]));
sidevars2D.push_back(new Variable2D(*sidevarsagain[3], *sidevarsagain[3]));
sidevars2D.push_back(new Variable2D(*sidevarsagain[7], *sidevarsagain[12]));
sidevars2D.push_back(new Variable2D(*sidevarsagain[8], *sidevarsagain[13]));
sidevars2D.push_back(new Variable2D(*sidevarsagain[9], *sidevarsagain[14]));
sidevars2D.push_back(new Variable2D(*sidevarsagain[10], *sidevarsagain[15]));
sidevars2D.push_back(new Variable2D(*sidevarsagain[11], *sidevarsagain[16]));
//vars2D.push_back(new Variable2D(vars[3],vars[0]));
if(doCCQENuValidation)
{
std::cerr << "Detected that tree name is CCQENu. Making validation histograms.\n";
vars.push_back(new Variable("pzmu", "p_{||, #mu} [GeV/c]", dansPzBins, &CVUniverse::GetMuonPz, &CVUniverse::GetMuonPzTrue));
vars.push_back(new Variable("Emu", "E_{#mu} [GeV]", robsEmuBins, &CVUniverse::GetEmuGeV, &CVUniverse::GetElepTrueGeV));
vars2D.push_back(new Variable2D(*vars[1], *vars[0]));
}
//vars2D.push_back(new Variable2D("q3_vs_pT", "q3_vs_pT;q3;pT", &CVUniverse::GetMuonPT, &CVUniverse::Getq3, dansPTBins,dansPTBins));
//This is where your list of Studies go for PerMichel variables -> Accessed through MichelEvent
std::vector<Study*> studies;
std::vector<Study*> sideband_studies;
std::function<double(const CVUniverse&, const MichelEvent&)> getq3 = [](const CVUniverse& univ, const MichelEvent& evt)
{
double q3 = evt.q3_reco;
//std::cout << "Printing q3 " << q3 << std::endl;
return q3;
};
std::function<double(const CVUniverse&, const MichelEvent&)> getpT = [](const CVUniverse& univ, const MichelEvent& evt)
{
double pT = evt.pT_reco;
//std::cout << "Printing pT " << pT << std::endl;
return pT;
};
std::function<double(const CVUniverse&, const MichelEvent&)> getEavail = [](const CVUniverse& univ, const MichelEvent& evt)
{
double eavail = evt.eavail_reco;
//std::cout << "Printing eavail " << eavail << std::endl;
return eavail;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> delta_t = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
int evttype = evt.eventtype;
double micheltime = evt.m_nmichels[whichMichel].time;
double vtxtime = univ.GetVertex().t();
double deltat = (micheltime - vtxtime/1000.); //hopefully this is in microseconds (mus)
//if (evttype == 1) return deltat;
//else return -9999.;
return deltat;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> permichel_range = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
//if (evt.eventtype == 1) return micheldist;
//else return -9999.;
return micheldist;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> pertruepimichel_range = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
if (evt.m_nmichels[whichMichel].true_parentpdg == 211) return micheldist;
else return -9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> permichel_tpi = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
if (evt.m_nmichels[whichMichel].true_parentpdg == 211) return evt.m_nmichels[whichMichel].pionKE;
else return 9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> overlay_vtx_tpi = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
double KE = evt.m_nmichels[whichMichel].pionKE;
int matchtype = evt.m_nmichels[whichMichel].BestMatch;
int overlayfrac = evt.m_nmichels[whichMichel].overlay_fraction;
if (overlayfrac > .5 && (matchtype == 1 || matchtype == 2 )) return KE;
else return -9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> overlay_clus_tpi = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
double KE = evt.m_nmichels[whichMichel].pionKE;
int matchtype = evt.m_nmichels[whichMichel].BestMatch;
int overlayfrac = evt.m_nmichels[whichMichel].overlay_fraction;
if (overlayfrac > .5 && (matchtype == 3 || matchtype == 4 )) return KE;
else return -9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> truemichel_goodvtx_tpi = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
int overlayfrac = evt.m_nmichels[whichMichel].overlay_fraction;
double KE = evt.m_nmichels[whichMichel].pionKE;
int matchtype = evt.m_nmichels[whichMichel].BestMatch;
int trueEnd = evt.m_nmichels[whichMichel].trueEndpoint;
int recoEnd = evt.m_nmichels[whichMichel].recoEndpoint;
if (overlayfrac < .5 && trueEnd == recoEnd && (matchtype == 1 || matchtype == 2)) return KE;
else return -9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> truemichel_goodclus_tpi = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
int overlayfrac = evt.m_nmichels[whichMichel].overlay_fraction;
double KE = evt.m_nmichels[whichMichel].pionKE;
int matchtype = evt.m_nmichels[whichMichel].BestMatch;
int trueEnd = evt.m_nmichels[whichMichel].trueEndpoint;
int recoEnd = evt.m_nmichels[whichMichel].recoEndpoint;
if (overlayfrac < .5 && trueEnd == recoEnd && (matchtype == 3 || matchtype == 4)) return KE;
else return -9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> truemichel_badvtx_tpi = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
double KE = evt.m_nmichels[whichMichel].pionKE;
int overlayfrac = evt.m_nmichels[whichMichel].overlay_fraction;
int matchtype = evt.m_nmichels[whichMichel].BestMatch;
int trueEnd = evt.m_nmichels[whichMichel].trueEndpoint;
int recoEnd = evt.m_nmichels[whichMichel].recoEndpoint;
if (overlayfrac < .5 && trueEnd != recoEnd && (matchtype == 1 || matchtype == 2))
{
//univ.PrintArachneLink();
//std::cout << "Printing Michel Time for bad VTX match type " << evt.m_nmichels[whichMichel].time << std::endl;
return KE;
}
else return -9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> truemichel_badclus_tpi = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
double KE = evt.m_nmichels[whichMichel].pionKE;
int overlayfrac = evt.m_nmichels[whichMichel].overlay_fraction;
int matchtype = evt.m_nmichels[whichMichel].BestMatch;
int trueEnd = evt.m_nmichels[whichMichel].trueEndpoint;
int recoEnd = evt.m_nmichels[whichMichel].recoEndpoint;
if (overlayfrac < .5 && trueEnd != recoEnd && (matchtype == 3 || matchtype == 4))
{
//univ.PrintArachneLink();
//std::cout << "Printing Michel Time for bad CLUS match type " << evt.m_nmichels[whichMichel].time << std::endl;
return KE;
}
else return -9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> overlay_vtx_range = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
int matchtype = evt.m_nmichels[whichMichel].BestMatch;
int overlayfrac = evt.m_nmichels[whichMichel].overlay_fraction;
if (overlayfrac > .5 && (matchtype == 1 || matchtype == 2 )) return micheldist;
else return -9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> overlay_clus_range = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
int matchtype = evt.m_nmichels[whichMichel].BestMatch;
int overlayfrac = evt.m_nmichels[whichMichel].overlay_fraction;
if (overlayfrac > .5 && (matchtype == 3 || matchtype == 4 )) return micheldist;
else return -9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> truemichel_goodvtx_range = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
int overlayfrac = evt.m_nmichels[whichMichel].overlay_fraction;
int matchtype = evt.m_nmichels[whichMichel].BestMatch;
int trueEnd = evt.m_nmichels[whichMichel].trueEndpoint;
int recoEnd = evt.m_nmichels[whichMichel].recoEndpoint;
if (overlayfrac < .5 && trueEnd == recoEnd && (matchtype == 1 || matchtype == 2)) return micheldist;
else return -9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> truemichel_goodclus_range = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
int overlayfrac = evt.m_nmichels[whichMichel].overlay_fraction;
int matchtype = evt.m_nmichels[whichMichel].BestMatch;
int trueEnd = evt.m_nmichels[whichMichel].trueEndpoint;
int recoEnd = evt.m_nmichels[whichMichel].recoEndpoint;
if (overlayfrac < .5 && trueEnd == recoEnd && (matchtype == 3 || matchtype == 4)) return micheldist;
else return -9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> truemichel_badvtx_range = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
int overlayfrac = evt.m_nmichels[whichMichel].overlay_fraction;
int matchtype = evt.m_nmichels[whichMichel].BestMatch;
int trueEnd = evt.m_nmichels[whichMichel].trueEndpoint;
int recoEnd = evt.m_nmichels[whichMichel].recoEndpoint;
if (overlayfrac < .5 && trueEnd != recoEnd && (matchtype == 1 || matchtype == 2))
{
// univ.PrintArachneLink();
// std::cout << "Printing Michel Time for bad VTX match type " << evt.m_nmichels[whichMichel].time << std::endl;
return micheldist;
}
else return -9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> truemichel_badclus_range = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
int overlayfrac = evt.m_nmichels[whichMichel].overlay_fraction;
int matchtype = evt.m_nmichels[whichMichel].BestMatch;
int trueEnd = evt.m_nmichels[whichMichel].trueEndpoint;
int recoEnd = evt.m_nmichels[whichMichel].recoEndpoint;
if (overlayfrac < .5 && trueEnd != recoEnd && (matchtype == 3 || matchtype == 4))
{
// univ.PrintArachneLink();
// std::cout << "Printing Michel Time for bad CLUS match type " << evt.m_nmichels[whichMichel].time << std::endl;
return micheldist;
}
else return -9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> michel_XZ = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double twoDdist = evt.m_nmichels[whichMichel].best_XZ;
return twoDdist;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> michel_UZ = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double twoDdist = evt.m_nmichels[whichMichel].best_UZ;
return twoDdist;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> michel_VZ = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double twoDdist = evt.m_nmichels[whichMichel].best_VZ;
return twoDdist;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> michel_XZ_upvtx = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double twoDdist = evt.m_nmichels[whichMichel].up_to_vertex_XZ;
int trueend = evt.m_nmichels[whichMichel].trueEndpoint;
if (trueend == 1) return twoDdist;
else return 9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> michel_UZ_upvtx = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double twoDdist = evt.m_nmichels[whichMichel].up_to_vertex_UZ;
int trueend = evt.m_nmichels[whichMichel].trueEndpoint;
if (trueend == 1) return twoDdist;
else return 9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> michel_VZ_upvtx = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double twoDdist = evt.m_nmichels[whichMichel].up_to_vertex_VZ;
int trueend = evt.m_nmichels[whichMichel].trueEndpoint;
if (trueend == 1) return twoDdist;
else return 9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> michel_XZ_upclus = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double twoDdist = evt.m_nmichels[whichMichel].up_to_clus_XZ;
int trueend = evt.m_nmichels[whichMichel].trueEndpoint;
if (trueend == 1) return twoDdist;
else return 9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> michel_UZ_upclus = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double twoDdist = evt.m_nmichels[whichMichel].up_to_clus_UZ;
int trueend = evt.m_nmichels[whichMichel].trueEndpoint;
if (trueend == 1) return twoDdist;
else return 9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> michel_VZ_upclus = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double twoDdist = evt.m_nmichels[whichMichel].up_to_clus_VZ;
int trueend = evt.m_nmichels[whichMichel].trueEndpoint;
if (trueend == 1) return twoDdist;
else return 9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> michel_XZ_downclus = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double twoDdist = evt.m_nmichels[whichMichel].down_to_clus_XZ;
int trueend = evt.m_nmichels[whichMichel].trueEndpoint;
if (trueend == 2) return twoDdist;
else return 9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> michel_UZ_downclus = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double twoDdist = evt.m_nmichels[whichMichel].down_to_clus_UZ;
int trueend = evt.m_nmichels[whichMichel].trueEndpoint;
if (trueend == 2) return twoDdist;
else return 9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> michel_VZ_downclus = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double twoDdist = evt.m_nmichels[whichMichel].down_to_clus_VZ;
int trueend = evt.m_nmichels[whichMichel].trueEndpoint;
if (trueend == 2) return twoDdist;
else return 9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> michel_XZ_downvtx = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double twoDdist = evt.m_nmichels[whichMichel].down_to_vertex_XZ;
int trueend = evt.m_nmichels[whichMichel].trueEndpoint;
if (trueend == 2) return twoDdist;
else return 9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> michel_UZ_downvtx = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double twoDdist = evt.m_nmichels[whichMichel].down_to_vertex_UZ;
int trueend = evt.m_nmichels[whichMichel].trueEndpoint;
if (trueend == 2) return twoDdist;
else return 9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> michel_VZ_downvtx = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double twoDdist = evt.m_nmichels[whichMichel].down_to_vertex_VZ;
int trueend = evt.m_nmichels[whichMichel].trueEndpoint;
if (trueend == 2) return twoDdist;
else return 9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> pion_angle = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double angle = evt.m_nmichels[whichMichel].best_angle;
return cos(angle);
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> pion_angle_range1 = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double angle = evt.m_nmichels[whichMichel].best_angle;
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
if (micheldist <= 150.) return cos(angle);
else return 9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> pion_angle_range2 = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{
double angle = evt.m_nmichels[whichMichel].best_angle;
double micheldist = evt.m_nmichels[whichMichel].Best3Ddist;
if (micheldist >150. && micheldist <= 250.) return cos(angle);
else return 9999.;
};
std::function<double(const CVUniverse&, const MichelEvent&, const int)> pion_angle_range3 = [](const CVUniverse& univ, const MichelEvent& evt, const int whichMichel)
{