HighLAND
CategoriesUtils.cxx
1 #include "ND280AnalysisUtils.hxx"
2 #include "BasicUtils.hxx"
3 #include "CategoriesUtils.hxx"
4 #include "TreeConverterUtils.hxx"
5 
6 
7 // PUT HERE ONLY COMMON CATEGORIES!
8 // YOU CAN CREATE YOUR OWN CATEGORY ON THE FLY IN A MACRO, LIKE THIS:
9 // draw.ChangeTrackCategory(<microtree_variable>, <ntypes>, <cat_types>, <cat_codes>, <cat_colors>);
10 // where the last 3 arguments are vectors similar to those defined here below and ntypes is their lenght
11 
12 // "no truth" and "sand #mu" types are added automatically to each category (see AddCategory)
13 
14 std::string nutype_types[] = {"#nu_{#mu}", "#bar{#nu}_{#mu}", "#nu_{e}", "#bar{#nu}_{e}", NAMEOTHER};
15 int nutype_codes[] = {14 , -14 , 12 , -12 , CATOTHER};
16 int nutype_colors[] = {2 , 3 , 4 , 5 , COLOTHER};
17 const int NNUTYPE = sizeof(nutype_types)/sizeof(nutype_types[0]);
18 
19 std::string part_types[] = {"#mu^{-}", "e^{-}", "#pi^{-}", "#mu^{+}", "e^{+}", "#pi^{+}", "p" , NAMEOTHER};
20 int pdg[] = {13 , 11 , -211 , -13 , -11 , 211 , 2212, CATOTHER};
21 int part_colors[] = {2 , 3 , 4 , 7 , 6 , 31 , 8 , COLOTHER};
22 const int NPART = sizeof(part_types)/sizeof(part_types[0]);
23 
24 std::string parent_types[] = {"#mu^{-}", "e^{-}", "#pi^{-}", "#mu^{+}", "e^{+}", "#pi^{+}", "p", "nu", "gamma", "#pi^{0}", "n" , NAMEOTHER};
25 int parent_pdg[] = {13 , 11 , -211 , -13 , -11 , 211 , 2212, 0 , 22 , 111 , 2112, CATOTHER};
26 int parent_colors[] = {2 , 3 , 4 , 7 , 6 , 31 , 1 , 8 , 9 , 40 , 95 , COLOTHER};
27 const int NPARENT = sizeof(parent_types)/sizeof(parent_types[0]);
28 
29 std::string reac_types[] = {"CCQE", NAME2P2H, "RES", "DIS", "COH", "NC", "#bar{#nu}_{#mu}", "#nu_{e}, #bar{#nu}_{e}", NAMEOTHER, NAMEOUTFV};
30 int reac_codes[] = {0 , CAT2P2H , 1 , 2 , 3 , 4 , 5 , 6 , CATOTHER , CATOUTFV};
31 int reac_colors[] = {2 , COL2P2H , 3 , 4 , 7 , 6 , 31 , 65 , COLOTHER , COLOUTFV};
32 const int NREAC = sizeof(reac_types)/sizeof(reac_types[0]);
33 
34 std::string reacCC_types[] = {"CC", "NC", "#bar{#nu}_{#mu}", "#nu_{e}, #bar{#nu}_{e}", NAMEOTHER, NAMEOUTFV};
35 int reacCC_codes[] = {1 , 4 , 5 , 6 , CATOTHER , CATOUTFV};
36 int reacCC_colors[] = {2 , 6 , 31 , 65 , COLOTHER , COLOUTFV};
37 const int NREACCC = sizeof(reacCC_types)/sizeof(reacCC_types[0]);
38 
39 std::string reacnofv_types[] = {"CCQE", NAME2P2H, "RES", "DIS", "COH", "NC", "#bar{#nu}_{#mu}", "#nu_{e}, #bar{#nu}_{e}", NAMEOTHER};
40 int reacnofv_codes[] = {0 , CAT2P2H , 1 , 2 , 3 , 4 , 5 , 6 , CATOTHER};
41 int reacnofv_colors[] = {2 , COL2P2H , 3 , 4 , 7 , 6 , 31 , 65 , COLOTHER};
42 const int NREACNOFV = sizeof(reacnofv_types)/sizeof(reacnofv_types[0]);
43 
44 std::string topology_types[] = {"CC-0#pi", "CC-1#pi^{+}", "CC-other", NAMEBKG, NAMEOUTFV};
45 int topology_codes[] = {0 , 1 , 2 , 3 , CATOUTFV};
46 int topology_colors[] = {2 , 4 , 7 , COLBKG , COLOUTFV};
47 const int NTOPOLOGY = sizeof(topology_types)/sizeof(topology_types[0]);
48 
49 std::string topology_no1pi_types[] = {"CC-0#pi", "CC-N#pi", NAMEBKG, NAMEOUTFV};
50 int topology_no1pi_codes[] = {0 , 1 , 2 , CATOUTFV};
51 int topology_no1pi_colors[] = {2 , 3 , COLBKG , COLOUTFV};
52 const int NTOPOLOGYNO1PI = sizeof(topology_no1pi_types)/sizeof(topology_no1pi_types[0]);
53 
54 std::string topology_withpi0_types[] = {"CC-0#pi", "CC-1#pi", "CC-X#pi^{0}", "CC-other", NAMEBKG, NAMEOUTFV};
55 int topology_withpi0_codes[] = {0 , 1 , 2 , 3 , 4 , CATOUTFV};
56 int topology_withpi0_colors[] = {2 , 4 , 31 , 65 , COLBKG , COLOUTFV};
57 const int NTOPOLOGYWITHPI0 = sizeof(topology_withpi0_types)/sizeof(topology_withpi0_types[0]);
58 
59 std::string topology_ccpizero_types[] = {"CC-1#pi^{0}", "CC-#pi^{0}+X", "CC-other", NAMEBKG, NAMEOUTFV};
60 int topology_ccpizero_codes[] = {0 , 1 , 2 , 3 , CATOUTFV};
61 int topology_ccpizero_colors[] = {2 , 4 , 7 , COLBKG , COLOUTFV};
62 const int NTOPOLOGYCCPIZERO = sizeof(topology_ccpizero_types)/sizeof(topology_ccpizero_types[0]);
63 
64 std::string mectopology_types[] = {"CC-0#pi-0p", "CC-0#pi-1p", "CC-0#pi-Np", "CC-1#pi^{+}", "CC-other", NAMEBKG, NAMEOUTFV};
65 int mectopology_codes[] = {0 , 1 , 2 , 3 , 4 , 5 , CATOUTFV};
66 int mectopology_colors[] = {2 , 3 , 4 , 7 , 31 , COLBKG , COLOUTFV};
67 const int NMECTOPOLOGY = sizeof(mectopology_types)/sizeof(mectopology_types[0]);
68 
69 std::string target_types[] = {"Carbon", "Oxygen", "Hydrogen", "Aluminium", "Iron", "Copper", "Zinc", "Lead", NAMEOTHER, "#splitline{coh on H}{bugzilla 1056}"};
70 int target_codes[] = {6 , 8 , 1 , 13 , 26 , 29 , 30 , 82 , CATOTHER, 0 };
71 int target_colors[] = {2 , 4 , 3 , 7 , 6 , 51 , 31 , 5 , COLOTHER, 40};
72 const int NTARGET = sizeof(target_types)/sizeof(target_types[0]);
73 
74 std::string sense_types[] = {"forward-correct", "forward-wrong", "backward-correct", "backward-wrong", "no z component", "unknown"};
75 int sense_codes[] = {0 , 1 , 2 , 3 , 4 , 5 };
76 int sense_colors[] = {2 , 3 , 4 , 7 , 6 , 31 };
77 
78 //TMP: keep the enumeration from the AnaTrueVertex::Detector (the original oaAnalysis one)
79 std::string detector_types[] = {"TPC1", "TPC2", "TPC3", "FGD1", "FGD2", "DsECAL", "BrECAL", "P0DECAL", "P0D", "SMRD", NAMEOTHER};
80 int detector_codes[] = {6, //SubDetId::kTPC1,
81  7, //SubDetId::kTPC2,
82  8, //SubDetId::kTPC3,
83  0, //SubDetId::kFGD1,
84  1, //SubDetId::kFGD2,
85  3, //SubDetId::kDSECAL,
86  4, //SubDetId::kTECAL,
87  5, //SubDetId::kPECAL,
88  2, //SubDetId::kP0D,
89  9, //SubDetId::kSMRD,
90  CATOTHER
91  };
92 int detector_colors[] = {2, 3, 4, 7, 6, 31, 51, 1, 65, 9, COLOTHER};
93 const int NDETECTOR = sizeof(detector_types)/sizeof(detector_types[0]);
94 
95 std::string momrelerr_types[] = {"0<momrelerr<0.05", "0.05<momrelerr<0.1", "0.1<momrelerr<0.2", "0.2<momrelerr<0.3", "0.3<momrelerr<0.4", "0.4<momrelerr<0.5", "0.5<momrelerr<0.6", "0.6<momrelerr<0.8", "0.8<momrelerr<1"};
96 int momrelerr_codes[] = { 0 , 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8};
97 int momrelerr_colors[] = { 1 , 2 , 4 , 3 , 6 , 7 , 31 , 5 , 8};
98 
99 
100 //********************************************************************
101 void anaUtils::AddStandardCategories(const std::string& prefix){
102 //********************************************************************
103 
104  // Add standard track categories for color drawing
105  _categ->AddCategory(prefix+"particle", NPART, part_types, pdg, part_colors);
106  _categ->AddCategory(prefix+"parent", NPARENT, parent_types, parent_pdg, parent_colors);
107  _categ->AddCategory(prefix+"gparent", NPARENT, parent_types, parent_pdg, parent_colors);
108  _categ->AddCategory(prefix+"primary", NPARENT, parent_types, parent_pdg, parent_colors);
109 
110  _categ->AddCategory(prefix+"nutype", NNUTYPE, nutype_types, nutype_codes, nutype_colors);
111  _categ->AddCategory(prefix+"nuparent", NPARENT, parent_types, parent_pdg, parent_colors);
112  _categ->AddCategory(prefix+"target", NTARGET, target_types, target_codes, target_colors);
113  _categ->AddCategory(prefix+"detector", NDETECTOR, detector_types, detector_codes, detector_colors);
114 
115  _categ->AddCategory(prefix+"reaction", NREAC, reac_types, reac_codes, reac_colors);
116  _categ->AddCategory(prefix+"reactionCC", NREACCC, reacCC_types, reacCC_codes, reacCC_colors);
117  _categ->AddCategory(prefix+"reactionnofv", NREACNOFV, reacnofv_types, reacnofv_codes, reacnofv_colors);
118  _categ->AddCategory(prefix+"topology", NTOPOLOGY, topology_types, topology_codes, topology_colors);
119  _categ->AddCategory(prefix+"topology_no1pi", NTOPOLOGYNO1PI, topology_no1pi_types, topology_no1pi_codes, topology_no1pi_colors);
120  _categ->AddCategory(prefix+"topology_withpi0", NTOPOLOGYWITHPI0, topology_withpi0_types, topology_withpi0_codes, topology_withpi0_colors);
121  _categ->AddCategory(prefix+"topology_ccpizero", NTOPOLOGYCCPIZERO, topology_ccpizero_types, topology_ccpizero_codes, topology_ccpizero_colors);
122  _categ->AddCategory(prefix+"mectopology", NMECTOPOLOGY, mectopology_types, mectopology_codes, mectopology_colors);
123 
124 }
125 
126 //********************************************************************
127 void anaUtils::AddStandardAntiNumuCategories(const std::string& prefix){
128 //********************************************************************
129 
130  // replace some naming
131  reac_types[6] = "#nu_{#mu}";
132  reacCC_types[2] = "#nu_{#mu}";
133  reacnofv_types[6] = "#nu_{#mu}";
134  mectopology_types[0] = "CC-0#pi-0n";
135  mectopology_types[1] = "CC-0#pi-1n";
136  mectopology_types[2] = "CC-0#pi-Nn";
137  mectopology_types[3] = "CC-1#pi^{-}";
138  topology_types[1] = "CC-1#pi^{-}";
139 
140  // Add standard anti-numu track categories for color drawing
141  AddStandardCategories(prefix);
142 }
143 
144 //**************************************************
145 Int_t anaUtils::GetTopology(const AnaTrueVertex& trueVertex, const SubDetId::SubDetEnum det, bool IsAntinu){
146 //**************************************************
147 
148  /* Classify reaction types
149  -1 = no true vertex
150  0 = CC 0pi = mu + X nucleons (+ 0 mesons)
151  1 = CC 1pi+(-) = mu + X nucleons + 1 pion +(-)
152  2 = CC other
153  3 = BKG: not (anti-)numu, NC
154  7 = out of FV
155  */
156 
157  // out of FGD1 FV
158  if ( ! anaUtils::InFiducialVolume(det,trueVertex.Position)) return CATOUTFV;
159 
160  /// primary particles from the true vertex
161  Int_t Nmuon = trueVertex.NPrimaryParticles[ParticleId::kMuon];
162  Int_t Nantimuon = trueVertex.NPrimaryParticles[ParticleId::kAntiMuon];
163  Int_t Nmeson = trueVertex.NPrimaryParticles[ParticleId::kMesons];
164  Int_t Npipos = trueVertex.NPrimaryParticles[ParticleId::kPiPos];
165  Int_t Npineg = trueVertex.NPrimaryParticles[ParticleId::kPiNeg];
166 
167  // non numu CC, i.e. BKG
168  if ( ! IsAntinu && Nmuon <= 0) return BKG; // BKG
169  if (IsAntinu && Nantimuon <= 0) return BKG; // BKG
170 
171  // CC 0pi, i.e. mu & 0 mesons
172  if (Nmeson == 0) return CC_0pi_0meson;
173 
174  // CC 1pi+, i.e. mu & 1 pi & no other mesons
175  if (Nmeson == 1) {
176  if ( ! IsAntinu && Npipos == 1) return CC_1pi_0meson;
177  else if (IsAntinu && Npineg == 1) return CC_1pi_0meson;
178  }
179 
180  // CC other
181  return CC_other;
182 }
183 
184 //**************************************************
185 Int_t anaUtils::GetTopology_no1pi(const AnaTrueVertex& trueVertex, const SubDetId::SubDetEnum det, bool IsAntinu){
186 //**************************************************
187 
188  /* Classify reaction types
189  -1 = no true vertex associated to track
190  0 = CC 0pi = mu + X nucleons (+ 0 mesons)
191  1 = CC other
192  2 = BKG: not (anti-)numu, NC
193  7 = out of FV
194  */
195 
196  // get GetTopology for numu
197  Int_t topo = GetTopology(trueVertex,det,IsAntinu);
198 
199  // merge CC_other and CC_1pi_0meson in the same category
200  if (topo == CC_other) return 1;
201  else if (topo == BKG) return 2;
202  else return topo;
203 }
204 
205 //**************************************************
206 Int_t anaUtils::GetTopology_withpi0(const AnaTrueVertex& trueVertex, const SubDetId::SubDetEnum det, bool IsAntinu){
207 //**************************************************
208 
209  /* Classify reaction types
210  -1 = no true vertex
211  0 = CC 0pi = mu + X nucleons (+ 0 mesons)
212  1 = CC 1pi+(-) = mu + X nucleons + 1 pion +(-)
213  2 = CC Npi0 = mu + Npi0 + X
214  3 = CC other
215  4 = BKG: not (anti-)numu, NC
216  7 = out of FV
217  */
218 
219  Int_t Npi0 = trueVertex.NPrimaryParticles[ParticleId::kPiZero];
220 
221  // get GetTopology for (anti-)numu
222  Int_t topo = GetTopology(trueVertex,det,IsAntinu);
223 
224  // split CC_other into with and w/o pi0
225  if (topo == CC_other) {
226  if (Npi0 > 0) return 2; // CC pi0 + X
227  else return 3;
228  }
229  else if (topo == BKG) return 4;
230  else return topo;
231 }
232 
233 //**************************************************
234 Int_t anaUtils::GetTopologyCCPiZero(const AnaTrueVertex& trueVertex, const SubDetId::SubDetEnum det, bool IsAntinu){
235 //**************************************************
236 
237  /* Classify the topology type for nu-mu CC pi-zero analysis
238  -1 = no true vertex
239  0 = CC 1pi0 = mu + X nucleons + 1 pi0
240  1 = CC Npi0 + X = mu + N>0 pi0 + X
241  2 = CC other
242  3 = BKG: not (anti-)numu CC
243  CATOUTFV = out of FV
244  */
245 
246  // out of FGD1 FV
247  if ( ! anaUtils::InFiducialVolume(det,trueVertex.Position)) return CATOUTFV;
248 
249  Int_t Nmuon = trueVertex.NPrimaryParticles[ParticleId::kMuon];
250  Int_t Nantimuon = trueVertex.NPrimaryParticles[ParticleId::kAntiMuon];
251  Int_t Nmeson = trueVertex.NPrimaryParticles[ParticleId::kMesons];
252  Int_t Npi0 = trueVertex.NPrimaryParticles[ParticleId::kPiZero];
253 
254  // non (anti-)numu CC, i.e. BKG
255  if ( ! IsAntinu && Nmuon <= 0) return 3; // BKG
256  if (IsAntinu && Nantimuon <= 0) return 3; // BKG
257 
258  // CC 1pi0
259  if (Npi0 == 1 && Nmeson == 1) return 0;
260 
261  // CC pi0 + X
262  else if (Npi0 > 0) return 1;
263 
264  // CC other
265  else return 2;
266 }
267 
268 //**************************************************
269 Int_t anaUtils::GetMECTopology(const AnaTrueVertex& trueVertex, const SubDetId::SubDetEnum det, bool IsAntinu){
270 //**************************************************
271 
272  /* Classify reaction types
273  -1 = no true vertex
274  0 = CC 0pi 0p(n) = mu + 0 p(n) (+ 0 mesons)
275  0 = CC 0pi 1p(n) = mu + 1 p(n) (+ 0 mesons)
276  0 = CC 0pi Np(n) = mu + N>1 p(n) (+ 0 mesons)
277  3 = CC 1pi+(-) = mu + X nucleons + 1 pion +(-)
278  4 = CC other : CC 1pi0, CC 1pi-, CCNpi+/-/0
279  5 = BKG: not (anti-)numu, NC
280  7 = out of FV
281  */
282 
283  Int_t Nproton = trueVertex.NPrimaryParticles[ParticleId::kProton];
284  Int_t Nneutron = trueVertex.NPrimaryParticles[ParticleId::kNeutron];
285 
286  // get GetTopology for (anti-)numu
287  Int_t topo = GetTopology(trueVertex,det,IsAntinu);
288 
289  if (topo == CC_0pi_0meson) {
290  Int_t Nnucleon;
291  if ( ! IsAntinu) Nnucleon = Nproton;
292  else Nnucleon = Nneutron;
293  if (Nnucleon == 0) return 0; // CC 0pi-0proton
294  if (Nnucleon == 1) return 1; // CC 0pi-1proton
295  if (Nnucleon > 1) return 2; // CC 0pi-Nproton
296  }
297  else if (topo == CC_1pi_0meson) return 3;
298  else if (topo == CC_other) return 4;
299  else if (topo == BKG) return 5;
300  else if (topo == CATOUTFV) return CATOUTFV;
301 
302  else std::cout << "Error in getting topology codes!" << std::endl;
303  return -999;
304 }
305 
306 //********************************************************************
307 void anaUtils::FillTruthTreeCategories(const AnaTrueVertex& trueVertex, const SubDetId::SubDetEnum det, bool IsAntinu){
308 //********************************************************************
309  FillCategories(&trueVertex,det,IsAntinu);
310 }
311 
312 //********************************************************************
313 void anaUtils::FillTruthTreeCategories(const AnaTrueVertex& trueVertex, const std::string& prefix, const SubDetId::SubDetEnum det, bool IsAntinu){
314 //********************************************************************
315  FillCategories(&trueVertex,prefix,det,IsAntinu);
316 }
317 
318 //********************************************************************
319 void anaUtils::FillCategories(AnaEventB* event, AnaTrack* track, const SubDetId::SubDetEnum det, bool IsAntinu, bool useCATSAND){
320 //********************************************************************
321  FillCategories(event, track,"",det,IsAntinu, useCATSAND);
322 }
323 
324 //********************************************************************
325 void anaUtils::SetCategoriesDefaultCode(const std::string& prefix, const int code){
326 //********************************************************************
327 
328  _categ->SetCode(prefix + "nutype" , code);
329  _categ->SetCode(prefix + "nuparent" , code);
330  _categ->SetCode(prefix + "detector" , code);
331  _categ->SetCode(prefix + "target" , code);
332 
333  _categ->SetCode(prefix + "reaction" , code);
334  _categ->SetCode(prefix + "reactionCC" , code);
335  _categ->SetCode(prefix + "reactionnofv" , code);
336  _categ->SetCode(prefix + "topology" , code);
337  _categ->SetCode(prefix + "topology_no1pi" , code);
338  _categ->SetCode(prefix + "topology_withpi0" , code);
339  _categ->SetCode(prefix + "topology_ccpizero" , code);
340  _categ->SetCode(prefix + "mectopology" , code);
341 }
342 
343 //********************************************************************
344 void anaUtils::FillCategories(AnaEventB* event, AnaTrack* track, const std::string& prefix, const SubDetId::SubDetEnum det, bool IsAntinu, bool useCATSAND){
345 //********************************************************************
346 
347  if ( ! track) return;
348 
349  // Categories are previously initialized with -999 in CategoryManager::ResetCurrentCategories()
350 
351  // Check if this is sandMC (by run number) even if there is no truth association (see bug 1237)
352  if (useCATSAND && event->GetIsSandMC()) {
353  SetCategoriesDefaultCode(prefix, CATSAND);
354  _categ->SetCode(prefix + "particle", CATSAND);
355  _categ->SetCode(prefix + "parent" , CATSAND);
356  _categ->SetCode(prefix + "gparent" , CATSAND);
357  _categ->SetCode(prefix + "primary" , CATSAND);
358  return;
359  }
360 
361  // First fill with "no truth" value
362  SetCategoriesDefaultCode(prefix, CATNOTRUTH);
363  _categ->SetCode(prefix + "particle", CATNOTRUTH);
364  _categ->SetCode(prefix + "parent" , CATNOTRUTH);
365  _categ->SetCode(prefix + "gparent" , CATNOTRUTH);
366  _categ->SetCode(prefix + "primary" , CATNOTRUTH);
367 
368  // ----- Particle ------------------------------
369  if (track->TrueObject) {
370 
371  if (track->GetTrueParticle()->PDG != 0) {
372  AnaTrueParticle* trueTrack = static_cast<AnaTrueParticle*>(track->TrueObject);
373  AnaTrueParticleB* primary = anaUtils::GetTrueParticleByID(*event, trueTrack->PrimaryID);
374 
375  if (trueTrack->TrueVertex) {
376  // Must be called before setting the particle category because by default it assumes the particle is the true lepton
377  // Overwritten below when there is a candidate
378  FillCategories(trueTrack->TrueVertex,prefix,det,IsAntinu);
379  }
380 
381  _categ->SetCode(prefix + "particle", trueTrack->PDG, CATOTHER);
382  _categ->SetCode(prefix + "parent", trueTrack->ParentPDG , CATOTHER);
383  _categ->SetCode(prefix + "gparent", trueTrack->GParentPDG, CATOTHER);
384 
385  if (primary) {
386  _categ->SetCode(prefix + "primary", primary->PDG, CATOTHER);
387  }
388 
389  }
390  }
391 }
392 
393 //********************************************************************
394 void anaUtils::FillCategories(const AnaTrueVertexB* trueVertexB, const std::string& prefix, const SubDetId::SubDetEnum det, bool IsAntinu, bool IsSand){
395 //********************************************************************
396 
397  if ( ! trueVertexB) return;
398 
399  // Categories are previously initialized with -999 in CategoryManager::ResetCurrentCategories()
400 
401  // if IsSand, assign CATSAND to categories
402  if (IsSand){
403  SetCategoriesDefaultCode(prefix, CATSAND);
404  return;
405  }
406 
407  SetCategoriesDefaultCode(prefix, CATNOTRUTH);
408 
409  const AnaTrueVertex* trueVertex = static_cast<const AnaTrueVertex*>(trueVertexB);
410 
411  // this only has effect when called initialy for a truth vertex (for TruthTree)
412  _categ->SetCode(prefix + "particle", trueVertex->LeptonPDG ,CATOTHER);
413 
414  _categ->SetCode(prefix + "nutype", trueVertex->NuPDG ,CATOTHER);
415  _categ->SetCode(prefix + "nuparent", trueVertex->NuParentPDG ,CATOTHER);
416 
417  int Detector_tmp;
419  _categ->SetCode(prefix + "detector", Detector_tmp ,CATOTHER);
420  _categ->SetCode(prefix + "target", GetTargetCode(trueVertex) ,CATOTHER);
421 
422  _categ->SetCode(prefix + "reaction", GetReaction(*trueVertex,det,IsAntinu) ,CATOTHER);
423  _categ->SetCode(prefix + "reactionCC", GetReactionCC(*trueVertex,det,IsAntinu) ,CATOTHER);
424  _categ->SetCode(prefix + "reactionnofv", GetReactionNoFgdFv(*trueVertex,IsAntinu) ,CATOTHER);
425  _categ->SetCode(prefix + "topology", GetTopology(*trueVertex,det,IsAntinu) ); // the default value has to be a category
426  _categ->SetCode(prefix + "topology_no1pi", GetTopology_no1pi(*trueVertex,det,IsAntinu) );
427  _categ->SetCode(prefix + "topology_withpi0", GetTopology_withpi0(*trueVertex,det,IsAntinu) );
428  _categ->SetCode(prefix + "topology_ccpizero", GetTopologyCCPiZero(*trueVertex, det,IsAntinu) );
429  _categ->SetCode(prefix + "mectopology", GetMECTopology(*trueVertex,det,IsAntinu) );
430 }
431 
432 //********************************************************************
433 void anaUtils::FillCategories(const AnaTrueVertexB* trueVertexB, const SubDetId::SubDetEnum det, bool IsAntinu, bool IsSand){
434 //********************************************************************
435 
436  FillCategories(trueVertexB,"",det,IsAntinu,IsSand);
437 }
438 
439 //**************************************************
440 Int_t anaUtils::GetReactionNoFgdFv(const AnaTrueVertex& trueVertex, bool IsAntinu){
441 //**************************************************
442  /* Classify reaction types
443  -1 = no true vertex
444  0 = CCQE
445  1 = CC RES
446  2 = CC DIS
447  3 = CC COH
448  4 = NC
449  5 = (anti-)numu
450  6 = nue/anti-nue
451  999 = other... nutau?
452  -1 = no truth
453  */
454 
455  Int_t reac = trueVertex.ReacCode;
456  Int_t nutype = trueVertex.NuPDG;
457 
458  // nue/anti-nue BKG
459  if (abs(nutype) == 12) return 6;
460 
461  // (anti-)nu BKG
462  else if (IsAntinu && nutype == 14) return 5; // nu BKG if antinu
463  else if ( ! IsAntinu && nutype == -14) return 5; // antinu BKG if nu
464 
465  // (anti-)numu
466  else if (abs(nutype) == 14) {
467  if (abs(reac) == 1) return 0;
468  if (abs(reac) == 2 ) return CAT2P2H;
469  // in NuWro prod 5 2p2h code is 70 (Neut does not have 70 at all)
470  if (abs(reac) == 70 ) return CAT2P2H;
471  if (abs(reac) >10 && abs(reac)<14) return 1;
472  if (abs(reac) >16 && abs(reac)<30) return 2;
473  if (abs(reac) ==16) return 3;
474  if (abs(reac) >30) return 4;
475  }
476 
477  return CATOTHER; // nu-tau??
478 }
479 
480 //**************************************************
481 Int_t anaUtils::GetReaction(const AnaTrueVertex& trueVertex, const SubDetId::SubDetEnum det, bool IsAntinu){
482 //**************************************************
483 
484  /* Classify reaction types
485  7 = out of FV
486  else as GetReactionNoFgdFv
487  */
488 
489  // out of FV
490  if ( ! InFiducialVolume(det, trueVertex.Position)) return CATOUTFV;
491 
492  return GetReactionNoFgdFv(trueVertex,IsAntinu);
493 }
494 
495 //**************************************************
496 Int_t anaUtils::GetReactionCC(const AnaTrueVertex& trueVertex, const SubDetId::SubDetEnum det, bool IsAntinu){
497 //**************************************************
498 
499  /* Classify reaction types
500  1 = CC (reaction<4 || 2P2H)
501  else = GetReaction
502  */
503 
504  Int_t reac = GetReaction(trueVertex,det,IsAntinu);
505  if (reac >= 0 && reac < 4) return 1;
506  else if (reac == CAT2P2H) return 1;
507  else return reac;
508 }
509 
510 //**************************************************
512 //**************************************************
513  Int_t code = -1;
514 
515  if (vtx) {
516  // Nucleus PDG codes are 10LZZZAAAI, and we want to extract the Z value to
517  // identify the element.
518  if (vtx->TargetPDG == 2212) // fix for bug in oaAnalysis, see bugzilla 1015
519  code = 1;
520  else
521  code = (vtx->TargetPDG / 10000) % 1000;
522 
523  }
524 
525  return code;
526 }
527 
528 //**************************************************
529 // antinu methods: if in .hxx "#ifndef CategoriesUtils_h" doesn't work
530 //**************************************************
531  Int_t anaUtils::GetTopology_antinu(const AnaTrueVertex& trueVertex, const SubDetId::SubDetEnum det) {return GetTopology(trueVertex,det,true);}
532  Int_t anaUtils::GetTopology_no1pi_antinu(const AnaTrueVertex& trueVertex, const SubDetId::SubDetEnum det) {return GetTopology_no1pi(trueVertex,det,true);}
533  Int_t anaUtils::GetTopology_withpi0_antinu(const AnaTrueVertex& trueVertex,const SubDetId::SubDetEnum det){return GetTopology_withpi0(trueVertex,det,true);}
534  Int_t anaUtils::GetTopologyCCPiZero_antinu(const AnaTrueVertex& trueVertex,const SubDetId::SubDetEnum det){return GetTopologyCCPiZero(trueVertex,det,true);}
535  Int_t anaUtils::GetMECTopology_antinu(const AnaTrueVertex& trueVertex, const SubDetId::SubDetEnum det) {return GetMECTopology(trueVertex,det,true);}
536  Int_t anaUtils::GetReaction_antinu(const AnaTrueVertex& trueVertex, const SubDetId::SubDetEnum det) {return GetReaction(trueVertex,det,true);}
537  Int_t anaUtils::GetReactionCC_antinu(const AnaTrueVertex& trueVertex, const SubDetId::SubDetEnum det) {return GetReactionCC(trueVertex,det,true);}
538  Int_t anaUtils::GetReactionNoFgdFv_antinu(const AnaTrueVertex& trueVertex) {return GetReactionNoFgdFv(trueVertex,true);}
539 //**************************************************
540 
541 
AnaTrueVertexB * TrueVertex
Pointer to the AnaTrueVertexB of the interaction that created this AnaTrueParticleB.
Int_t GetTopology_withpi0(const AnaTrueVertex &trueVertex, const SubDetId::SubDetEnum det=SubDetId::kFGD1, bool IsAntinu=false)
Classify the topology type for nu-mu CC pi-zero analysis.
Representation of a global track.
Int_t NuPDG
The PDG code of the incoming neutrino.
Int_t LeptonPDG
The PDG code of the primary outgoing electron/muon.
void AddCategory(const std::string &name, int ntypes, std::string *names, int *codes, int *colors, bool multi=false, bool noWarning=false, bool addNOTRUTH=true, bool addSAND=true)
void SetCode(const std::string &categ, int code, int defaultcode=0)
Int_t NuParentPDG
Neutrino parent PDG code.
Definition: DataClasses.hxx:88
bool GetIsSandMC() const
Return whether this event is from Sand Monte Carlo or not.
Representation of a true Monte Carlo vertex.
Int_t GParentPDG
The PDG code of this particle&#39;s grandparent, or 0 if there is no grandparent.
AnaTrueObjectC * TrueObject
The link to the true oject that most likely generated this reconstructed object.
Int_t GetMECTopology(const AnaTrueVertex &trueVertex, const SubDetId::SubDetEnum det=SubDetId::kFGD1, bool IsAntinu=false)
Classify reaction topologies in special attention to MEC process.
Int_t TargetPDG
The PDG code of the target nucleus.
Definition: DataClasses.hxx:85
Int_t GetTopology(const AnaTrueVertex &trueVertex, const SubDetId::SubDetEnum det=SubDetId::kFGD1, bool IsAntinu=false)
Classify reaction topologies.
Int_t GetReaction(const AnaTrueVertex &trueVertex, const SubDetId::SubDetEnum det=SubDetId::kFGD1, bool IsAntinu=false)
Classify reaction types.
Int_t GetTopology_antinu(const AnaTrueVertex &trueVertex, const SubDetId::SubDetEnum det=SubDetId::kFGD1)
Classify reaction types for antinu.
Representation of a true Monte Carlo vertex.
Definition: DataClasses.hxx:50
Representation of a true Monte Carlo trajectory/particle.
void AddStandardAntiNumuCategories(const std::string &prefix="")
Add the standard categories only, given a prefix for their name.
Float_t Position[4]
The position the true interaction happened at.
SubDetEnum
Enumeration of all detector systems and subdetectors.
Definition: SubDetId.hxx:25
Int_t PDG
The PDG code of this particle.
unsigned long Detector
Int_t ParentPDG
The PDG code of this particle&#39;s immediate parent, or 0 if there is no parent.
void FillCategories(AnaEventB *event, AnaTrack *track, const std::string &prefix, const SubDetId::SubDetEnum det=SubDetId::kFGD1, bool IsAntinu=false, bool useCATSAND=true)
Fill the track categories for color drawing.
AnaTrueParticleB * GetTrueParticleByID(const AnaEventB &event, int ID)
Get the AnaTrueParticleB in the current spill with the given ID. Return NULL if it can&#39;t be found...
void AddStandardCategories(const std::string &prefix="")
Add the standard categories only, given a prefix for their name.
AnaTrueParticleB * GetTrueParticle() const
Return a casted version of the AnaTrueObject associated.
bool InFiducialVolume(SubDetId::SubDetEnum det, const Float_t *pos, const Float_t *FVdefmin, const Float_t *FVdefmax)
Representation of a true Monte Carlo trajectory/particle.
void ConvertBitFieldToTrueParticleDetEnum(unsigned long det, int &trueDet)
Get highland style numbering for true objects.
Int_t GetTargetCode(const AnaTrueVertex *trueVertex)
Get the code for filling the target PDG category.
Int_t NPrimaryParticles[Int_t(ParticleId::kLast)+1]
Array to count the outgoing primary particles of each type (.