HighLAND
ND280AnalysisUtils.cxx
1 #include "ND280AnalysisUtils.hxx"
2 #include <TVector3.h>
3 #include <stdio.h>
4 #include <math.h>
5 #include <iostream>
6 #include <typeinfo>
7 #include "BaseDataClasses.hxx"
8 #include "DetectorDefinition.hxx"
9 #include <ND280BeamBunching.hxx>
10 #include "PionInteractionSystematic.hxx"
11 #include "Parameters.hxx"
12 
13 ND280BeamBunching bunching;
14 
15 
16 //**************************************************
17 Int_t anaUtils::GetRunPeriod(Int_t run, Int_t subrun){
18 //**************************************************
19 
20  if ( run<6000) return 0;//run1 water
21  else if (run>=6462 && run<7664) return 1;//run2 water
22  else if (run>=7665 && run<7754) return 2;//run2 air
23  else if (run>=8000 && run<8550) return 3;//run3b air => mc associated to be the one of run2
24  else if (run>=8550 && run<8800) return 4;//run3c air
25  else if (run>=8983 && run<9414) return 5;//run4 water
26  else if (run>=9426 && run<9800) return 6;//run4 air
27  else if ((run>=10252 && run<10334) // + 10252_10 to 10333
28  || (run>=10519 && run<10522)) return 7;//run5 water + 10518_28 to 10521_13
29  else if (run>=10334 && run<10518) return 8;//run5 water (ANTINU) up to 10334_11 to 10518_9
30  else if (run == 10518){
31  if(subrun < 10) return 8; //run5 water (ANTINU)
32  else return 7; //run5 water
33  }
34  else if (run>=10828 && run<10954) return 13; //run6 air
35  else if (run == 10954){
36  if(subrun > -1 && subrun < 10) return 13; //run6 air
37  else return 9; // run6b air (ANTINU)
38  }
39  else if (run>=10955 && run<11240) return 9;//run6b air (ANTINU)
40  else if (run>=11305 && run<11443) return 13; //run6 air
41  else if (run>=11449 && run<11492) return 10;//run6c air (ANTINU)
42  else if (run>=11492 && run<11564) return 11;//run6d air (ANTINU)
43  else if (run>=11619 && run<11687) return 12;//run6e air (ANTINU)
44  else if (run == 11687){
45  if(subrun < 7) return 12;//run6e air (ANTINU)
46  else return 13; //run6 air
47  }
48 
49  else if (run>=90110000 && run<=90119999) return 0;//run1 water
50  else if (run>=90210000 && run<=90219999) return 1;//run2 water
51  else if (run>=90200000 && run<=90209999) return 2;//run2 air
52  else if (run>=90300000 && run<=90300015) return 3;//run3b air separated based on the pot of data (run3b/(run3b+run3c)=0.13542
53  else if (run>=90300016 && run<=90300110) return 4;//run3c air separated based on the pot of data
54  else if (run>=90410000 && run<=90419999) return 5;//run4 water
55  else if (run>=90400000 && run<=90409999) return 6;//run4 air
56  else if (run>=80510000 && run<=80519999) return 8;//run5 antinu-water
57  else if (run>=80600000 && run<=80600159) return 9;//run6b antinu-air
58  else if (run>=80600160 && run<=80600219) return 10;//run6c antinu-air - have to split Run 6 due to different flux tunings for the different parts
59  else if (run>=80600220 && run<=80600299) return 11;//run6d antinu-air
60  else if (run>=80600300 && run<=80600399) return 12;//run6e antinu-air
61  else if (run>=80600400 && run<=80609999) return 12;//run6e antinu-air - set this as default, to catch any files that are not currently at TRIUMF
62 
63  else if (run>=80307000 && run<=80307999) return 8;//sand muons antinu
64  else if (run>=90307000 && run<=90307999) return 4;//sand muons neut
65  else if (run>=92307000 && run<=92307999) return 4;//sand muons old neut
66 
67  // genie
68  else if (run>=91110000 && run<=91119999) return 0;//run1 water
69  else if (run>=91210000 && run<=91219999) return 1;//run2 water
70  else if (run>=91200000 && run<=91209999) return 2;//run2 air
71  else if (run>=91300000 && run<=91300015) return 3;//run3 air separated based on the pot of data (run3b/(run3b+run3c)=0.13542
72  else if (run>=91300016 && run<=91300110) return 4;//run3 air separated based on the pot of data (run3b/(run3b+run3c)=0.13542
73  else if (run>=91410000 && run<=91419999) return 5;//run4 water
74  else if (run>=91400000 && run<=91409999) return 6;//run4 air
75  else if (run>=81510000 && run<=81519999) return 8;//run5 antinu-water
76  else if (run>=81600000 && run<=81600159) return 9;//run6b antinu-air
77  else if (run>=81600160 && run<=81600219) return 10;//run6c antinu-air - have to split Run 6 due to different flux tunings for the different parts
78  else if (run>=81600220 && run<=81600299) return 11;//run6d antinu-air
79  else if (run>=81600300 && run<=81600399) return 12;//run6e antinu-air
80  else if (run>=81600400 && run<=81609999) return 12;//run6e antinu-air - set this as default, to catch any files that are not currently at TRIUMF
81  else return -1;
82 
83 }
84 
85 
86 //**************************************************
87 int anaUtils::GetSandMode(int run){
88 //**************************************************
89 
90  if (run>=80307000 && run<=80307999) return 0; // antineutrino flux
91  else if (run>=90307000 && run<=90307999) return 1; // neutrino flux
92  else if (run>=92307000 && run<=92307999) return 2; // neutrino flux, old NEUT
93 
94  else return -1;
95 }
96 
97 //*****************************************************************************
98 int anaUtils::GetFGDMichelElectrons(const AnaEventB& event, const SubDetId::SubDetEnum det, AnaFgdTimeBinB** arr, bool prod5Cut){
99 //*****************************************************************************
100 
101  int nSel = 0;
102 
103  Int_t ibunch = event.Bunch;
104  Float_t selTime = bunching.GetBunchCentralTime(event, ibunch);
105 
106  Float_t numWidthsReqInSorting = 4;
107  Float_t binBunchWidth = 25.;
108  Float_t binTimeOffsetToBunch = -13.;
109  Float_t minRecBinQ = 0.;
110  Int_t minHitsFGD1 = 6;
111  Int_t minHitsFGD2 = 5;
112 
113  if(prod5Cut){
114  minRecBinQ = 200.0;
115  minHitsFGD1 = 0;
116  minHitsFGD2 = 0;
117  }
118 
119  int index=-1;
120 
121  selTime += numWidthsReqInSorting*binBunchWidth;
122 
123  for( Int_t i = 0; i < event.nFgdTimeBins; i++ ){
124 
125  AnaFgdTimeBinB *FgdTimeBin = event.FgdTimeBins[i];
126 
127  if( det == SubDetId::kFGD1 && FgdTimeBin->NHits[0] == 0 ) continue;
128  if( det == SubDetId::kFGD2 && FgdTimeBin->NHits[1] == 0 ) continue;
129  if( det == SubDetId::kFGD && FgdTimeBin->NHits[0] == 0 && FgdTimeBin->NHits[1] == 0 ) continue;
130 
131  Float_t binTime = FgdTimeBin->MinTime - binTimeOffsetToBunch;
132 
133  bool inBunch = false;
134 
135  for( UInt_t k = 0 ; k < bunching.GetNBunches(event); k++ ){
136  Float_t time = bunching.GetBunchCentralTime(event, k);
137  if( TMath::Abs( binTime - time ) < numWidthsReqInSorting*binBunchWidth ) inBunch = true;
138  }
139 
140  if( inBunch ) continue;
141 
142  //find bin with min raw charge and greater than the required number of hits
143  if( binTime >= selTime && FgdTimeBin->RawChargeSum[0] > minRecBinQ && FgdTimeBin->NHits[0] > minHitsFGD1 && (det == SubDetId::kFGD || det == SubDetId::kFGD1) ){
144  index = i;
145  }
146  if( binTime >= selTime && FgdTimeBin->RawChargeSum[1] > minRecBinQ && FgdTimeBin->NHits[1] > minHitsFGD2 && (det == SubDetId::kFGD || det == SubDetId::kFGD2) ){
147  index = i;
148  }
149 
150  }
151  if(index > -1){
152  arr[nSel] = event.FgdTimeBins[index];
153  nSel++;
154  }
155 
156  return nSel;
157 }
158 
159 
160 // //*****************************************************************************
161 // bool anaUtils::IsFGD1Analysis(const AnaEventB& event){
162 // //*****************************************************************************
163 
164 // if(!event.Summary) return false;
165 // if(!event.Summary->EventSample) return false;
166 
167 // if((event.Summary->EventSample >= SampleId::kFGD1NuMuCC
168 // && event.Summary->EventSample<=SampleId::kFGD1NuMuCCOther) ||
169 // (event.Summary->EventSample >= SampleId::kFGD1AntiNuMuCC
170 // && event.Summary->EventSample<=SampleId::kFGD1AntiNuMuCCnQE)||
171 // (event.Summary->EventSample >= SampleId::kFGD1NuMuBkgInAntiNuModeCC
172 // && event.Summary->EventSample<=SampleId::kFGD1NuMuBkgInAntiNuModeCCnQE))
173 // return true;
174 // return false;
175 // }
176 // //*****************************************************************************
177 // bool anaUtils::IsFGD2Analysis(const AnaEventB& event){
178 // //*****************************************************************************
179 // if(!event.Summary) return false;
180 // if(!event.Summary->EventSample) return false;
181 
182 // if((event.Summary->EventSample >= SampleId::kFGD2NuMuCC
183 // && event.Summary->EventSample<= SampleId::kFGD2NuMuCCOther) ||
184 // (event.Summary->EventSample >= SampleId::kFGD2AntiNuMuCC
185 // && event.Summary->EventSample<= SampleId::kFGD2AntiNuMuCCNTracks) ||
186 // (event.Summary->EventSample >= SampleId::kFGD2NuMuBkgInAntiNuModeCC
187 // && event.Summary->EventSample <= SampleId::kFGD2NuMuBkgInAntiNuModeCCnQE))
188 // return true;
189 // return false;
190 // }
191 
192 //*********************************************************
193 double anaUtils::GetNucleonsPerNucleus(massComponentEnum massComponent, int target_nucleons) {
194 //*********************************************************
195  // target_nucleons: 1 for protons, 2 for neutrons, 3 for both
196 
197  double nNucleonsPerNucleus = 0., nProtonsPerNucleus=0;
198 
199  // Areal densities of each element in FGDs (=>relative abundance), taken from TN 91 and TN 198 (weighted averaged for FGD2)
200  // those commented are the values used in TN 117 and TN 210, where 2nd and 3rd column are rounded fractions, i.e. divided by the total areal density
201  // element[4] = {Nprotons,
202  // Nnucleons (initialized later),
203  // 'typical' areal densities of XY modules in g/cm^2,
204  // weighted averaged areal densities of water modules in mg/cm^2}
205  double carbon[4] = {6, 0, 1.849, 421.0704607 }; // {6, 0, 0.8613, 0.1505};
206  double oxygen[4] = {8, 0, 0.0794, (12397./6.) }; // {8, 0, 0.0370, 0.7383};
207  double hydrogen[4] = {1, 0, 0.1579, (1759./6.) }; // {1, 0, 0.0736, 0.1048};
208  double titanium[4] = {22, 0, 0.0355, 0. }; // {22, 0, 0.0165, 0. };
209  double silicon[4] = {14, 0, 0.0218, 11 }; // {14, 0, 0.0101, 0.0039};
210  double nitrogen[4] = {7, 0, 0.0031, 0. }; // {7, 0, 0.0014, 0. };
211  double magnesium[4] = {12, 0, 0., 7. }; // {12, 0, 0., 0.0025};
212  // totArealDensityXY = 2.1467 g/cm^2 summing, while in TN 91 is reported 2.147 and in NIM FGDs 2.1463
213  double totArealDensityXY = carbon[2] + oxygen[2] + hydrogen[2] + titanium[2] + silicon[2] + nitrogen[2] + magnesium[2];
214  double totArealDensityWater = carbon[3] + oxygen[3] + hydrogen[3] + titanium[3] + silicon[3] + nitrogen[3] + magnesium[3];
215 
216  // Nnucleons with isotopic abundance
217  // values taken from www.ciaaw.org/isotopic-abundances.htm: Elements published in Pure and Applied Chemistry in 2011
218  // those commented are the values used in TN 117 and TN 210
219  carbon[1] = 12 * 0.9893 + 13 * 0.0107; // 12*0.989 + 13*0.011;
220  oxygen[1] = 16 * 0.99757 + 17 * 0.00038 + 18 * 0.00205; // 16*0.99762 + 17*0.00038 + 18*0.002;
221  hydrogen[1] = 1 * 0.999885 + 2 * 0.000115; // 1*0.99985 + 2*0.00015;
222  titanium[1] = 46 * 0.0825 + 47 * 0.0744 + 48 * 0.7372 + 49 * 0.0541 + 50 * 0.0518; // 46*0.08 + 47*0.075 + 48*0.738 + 49*0.055 + 50*0.054;
223  silicon[1] = 28 * 0.92223 + 29 * 0.04685 + 30 * 0.03092; // 28*0.9222 + 29*0.0468 + 30*0.0309;
224  nitrogen[1] = 14 * 0.996346 + 15 * 0.00364; // 14*0.99634 + 15*0.00366;
225  magnesium[1] = 24 * 0.7899 + 25 * 0.1000 + 26 * 0.1101; // 24*0.7899 + 25*0.1 + 26*0.11;
226 
227  if(massComponent == kFGD1 || massComponent == kFGD2xymodules || massComponent == kFGD2xylike) {
228 
229  // average nucleons/protons per nucleus in FGD1
230  nProtonsPerNucleus = carbon[0] * carbon[2] +
231  oxygen[0] * oxygen[2] +
232  hydrogen[0] * hydrogen[2] +
233  titanium[0] * titanium[2] +
234  silicon[0] * silicon[2] +
235  nitrogen[0] * nitrogen[2] +
236  magnesium[0] * magnesium[2] ;
237  nProtonsPerNucleus = nProtonsPerNucleus / totArealDensityXY;
238 
239  nNucleonsPerNucleus = carbon[1] * carbon[2] +
240  oxygen[1] * oxygen[2] +
241  hydrogen[1] * hydrogen[2] +
242  titanium[1] * titanium[2] +
243  silicon[1] * silicon[2] +
244  nitrogen[1] * nitrogen[2] +
245  magnesium[1] * magnesium[2] ;
246  nNucleonsPerNucleus = nNucleonsPerNucleus / totArealDensityXY;
247  }
248  else if(massComponent == kFGD2watermodules || massComponent == kFGD2waterlike) {
249 
250  // average nucleons/protons per nucleus in FGD2
251  nProtonsPerNucleus = carbon[0] * carbon[3] +
252  oxygen[0] * oxygen[3] +
253  hydrogen[0] * hydrogen[3] +
254  titanium[0] * titanium[3] +
255  silicon[0] * silicon[3] +
256  nitrogen[0] * nitrogen[3] +
257  magnesium[0] * magnesium[3] ;
258  nProtonsPerNucleus = nProtonsPerNucleus / totArealDensityWater;
259 
260  nNucleonsPerNucleus = carbon[1] * carbon[3] +
261  oxygen[1] * oxygen[3] +
262  hydrogen[1] * hydrogen[3] +
263  titanium[1] * titanium[3] +
264  silicon[1] * silicon[3] +
265  nitrogen[1] * nitrogen[3] +
266  magnesium[1] * magnesium[3] ;
267  nNucleonsPerNucleus = nNucleonsPerNucleus / totArealDensityWater;
268  }
269  else if (massComponent == kP0Dwater) {
270 
271  nProtonsPerNucleus = (8*1/3. + // Oxygen
272  1*2/3.); // Hydrogen
273 
274  nNucleonsPerNucleus = ((16*0.99762+17*0.00038+18*0.002)*1/3. + // Oxygen
275  (1*0.99985+2*0.00015)*2/3.); // Hydrogen
276  }
277  else {
278  std::cout << "Unknown target section selected." << std::endl;
279  return 0;
280  }
281 
282  if (target_nucleons == 1) return nProtonsPerNucleus;
283  else if(target_nucleons == 2) return (nNucleonsPerNucleus - nProtonsPerNucleus);
284  else if(target_nucleons == 3) return nNucleonsPerNucleus;
285  else {
286  std::cout << "Unknown target nucleons selected." << std::endl;
287  return 0;
288  }
289 }
290 
291 //*********************************************************
292 double anaUtils::GetAreaFV(massComponentEnum massComponent) {
293 //*********************************************************
294  TVector3 fgd1min = anaUtils::ArrayToTVector3(DetDef::fgd1min);
295  TVector3 fgd2min = anaUtils::ArrayToTVector3(DetDef::fgd2min);
296  TVector3 fgd1max = anaUtils::ArrayToTVector3(DetDef::fgd1max);
297  TVector3 fgd2max = anaUtils::ArrayToTVector3(DetDef::fgd2max);
298  TVector3 FVdefminFGD1 = anaUtils::ArrayToTVector3(FVDef::FVdefminFGD1);
299  TVector3 FVdefminFGD2 = anaUtils::ArrayToTVector3(FVDef::FVdefminFGD2);
300  TVector3 FVdefmaxFGD1 = anaUtils::ArrayToTVector3(FVDef::FVdefmaxFGD1);
301  TVector3 FVdefmaxFGD2 = anaUtils::ArrayToTVector3(FVDef::FVdefmaxFGD2);
302 
303  TVector3 FVfgd1 = (fgd1max - FVdefmaxFGD1)-(fgd1min + FVdefminFGD1);
304  TVector3 FVfgd2 = (fgd2max - FVdefmaxFGD2)-(fgd2min + FVdefminFGD2);
305 
306  double areaFV = 0;
307  if (massComponent == kFGD1) {
308  areaFV = FVfgd1.x() * FVfgd1.y() / 100.; // convert from mm^2 to cm^2
309  } else if (massComponent == kFGD2xymodules || massComponent == kFGD2watermodules || massComponent == kFGD2waterlike || massComponent == kFGD2xylike) {
310  areaFV = FVfgd2.x() * FVfgd2.y() / 100.; // convert from mm^2 to cm^2
311  }
312  else {
313  std::cout << "Unknown target section selected." << std::endl;
314  return 0;
315  }
316 
317  return areaFV;
318 }
319 
320 //*********************************************************
321 double anaUtils::GetNTargets(massComponentEnum massComponent, int target_nucleons) {
322 //*********************************************************
323  // target_nucleons: 1 for protons, 2 for neutrons, 3 for both
324 
325  // XY modules
326  double nXYModulesFGD1FV = 15.; // whole volume
327  double nXYModulesFGD2FV = 7.; // whole volume
328  // avoid dealing with air between modules by looking only at FVdef and taking the integer of the division
329  double discardedModulesFGD1 = (FVDef::FVdefminFGD1[2] + FVDef::FVdefmaxFGD1[2]) / DetDef::fgdXYModuleWidth;
330  double discardedModulesFGD2 = (FVDef::FVdefminFGD2[2] + FVDef::FVdefmaxFGD2[2]) / DetDef::fgdXYModuleWidth;
331  // multiply (and then divide) by 2 to allow approximation to half module
332  int nFGD1 = (int)(discardedModulesFGD1 * 2) ;
333  int nFGD2 = (int)(discardedModulesFGD2 * 2) ;
334  nXYModulesFGD1FV -= (double)(nFGD1/2.);
335  nXYModulesFGD2FV -= (double)(nFGD2/2.);
336 
337  // Water modules
338  double xyplusair = DetDef::fgdXYModuleWidth + DetDef::fgdXYAirWidth;
339  if (FVDef::FVdefminFGD2[2] > xyplusair || FVDef::FVdefmaxFGD2[2] > xyplusair)
340  std::cout << "ERROR: FGD2 FV doesn't include the 6 water modules --> GetNTargets will be wrong for FGD2!" << std::endl;
341  double nWaterModules = 6.;
342 
343  double massFV = 0.;
344  if(massComponent == kFGD1 || massComponent == kFGD2xymodules) {
345 
346  double nXYModulesFV = 0;
347  if (massComponent == kFGD1) {
348  nXYModulesFV = nXYModulesFGD1FV;
349  std::cout << "FGD1 (" << nXYModulesFV << " modules):" << std::endl;
350  } else if (massComponent == kFGD2xymodules) {
351  nXYModulesFV = nXYModulesFGD2FV;
352  std::cout << "FGD2 XY (" << nXYModulesFV << " modules):" << std::endl;
353  }
354 
355  // areal density from tab 1 in TN 91, even though the sum seems wrong (vertically 2.1463, horizontally 2.1467)
356  // but air is missing (which is 0.000258 from TN 122)
357  massFV = 2.147 * GetAreaFV(massComponent) * nXYModulesFV;
358  double density = 2.147 / (((double)DetDef::fgdXYModuleWidth + (double)DetDef::fgdXYAirWidth) / 10.); // convert mm to cm
359  std::cout << " density (included air): " << density << " g/cm^3" << std::endl;
360  }
361 
362  else if(massComponent == kFGD2watermodules) {
363  // areal density from tab 9 in TN 198 (weighted average of 'total'), as built
364  // in tab 5 of TN 122 (bottom three rows) it is instead 2.791, but it is the MC one
365  // (in both air is missing but from TN 198 it is only 0.000266011)
366  massFV = 2.79867 * GetAreaFV(massComponent) * nWaterModules;
367  std::cout << "FGD2 water modules (" << nWaterModules << " modules):" << std::endl;
368  double density = 2.79867 / (((double)DetDef::fgdWaterModuleWidth + (double)DetDef::fgdWaterAirWidth) / 10.); // convert mm to cm
369  std::cout << " density (included air): " << density << " g/cm^3" << std::endl;
370  }
371 
372  else if(massComponent == kFGD2waterlike) {
373  // areal density from tab 78 in TN 212, already for 6 modules, which is the MC one
374  // it corresponds to the as built one, from tab 9 and 10 in TN 198 (weighted average of 'water' + 'virtual water')
375  massFV = 13.838 * GetAreaFV(massComponent);
376  std::cout << "FGD2 water-like (" << nWaterModules << " modules):" << std::endl;
377  }
378 
379  else if(massComponent == kFGD2xylike) { // this is XY + virtual XY
380  // areal density from tab 78 in TN 212, already for 6 modules, which is the MC one
381  // it corresponds to the as built one, from tab 9 and 10 in TN 198 (weighted average of 'water' + 'virtual water')
382  massFV = 16.887 * GetAreaFV(massComponent);
383  std::cout << "FGD2 XY-like (mixed modules):" << std::endl;
384  }
385 
386  else if(massComponent == kFGD2) {
387  double Nwater = GetNTargets(kFGD2watermodules,target_nucleons);
388  double Nscint = GetNTargets(kFGD2xymodules,target_nucleons);
389  return (Nwater + Nscint);
390  }
391 
392  else if(massComponent == kFGDs) {
393  double Nfgd1 = GetNTargets(kFGD1,target_nucleons);
394  double Nfgd2 = GetNTargets(kFGD2,target_nucleons);
395  return (Nfgd1 + Nfgd2);
396  }
397 
398  else if (massComponent == kP0Dwater) {
399  // Water mass from TN82v2 and p0dasbuilt spreadsheet
400  massFV = 1902.*1.e3; // g
401  std::cout << "P0D:" << std::endl;
402  }
403 
404  else {
405  std::cout << "Unknown target section selected." << std::endl;
406  return 0;
407  }
408 
409  std::cout << " mass FV: " << massFV / 1000. << " Kg" << std::endl; // convert g to Kg
410 
411  // number of nucleons (assuming 1g/mol/nucleon)
412  double nNucleons = massFV * units::Avogadro;
413  if (target_nucleons == 3) return nNucleons;
414  double nNucleonsPerNucleus = GetNucleonsPerNucleus(massComponent, 3);
415  return nNucleons * GetNucleonsPerNucleus(massComponent, target_nucleons) / nNucleonsPerNucleus;
416 }
417 
418 //********************************************************************
419 void anaUtils::IncrementPOTBySpill(const AnaSpillB& spill, Header& header) {
420 //********************************************************************
421 
422  const AnaBeamB& beam = *spill.Beam;
423  const AnaEventInfoB& info = *spill.EventInfo;
424 
425  if (beam.POTSincePreviousSavedSpill > 1e+16) {
426  std::cout << "WARNING: POT is suspiciously large for run " << info.Run << ", subrun " << info.SubRun << ", event " << info.Event << ": " << beam.POTSincePreviousSavedSpill << ". POT will not be counted for this event." << std::endl;
427  return;
428  }
429 
430  if (beam.POTSincePreviousSavedSpill < 0) {
431  std::cout << "WARNING: POT is negative for run " << info.Run << ", subrun " << info.SubRun << ", event " << info.Event << ": " << beam.POTSincePreviousSavedSpill << ". POT will not be counted for this event." << std::endl;
432  return;
433  }
434 
435  // For real data
436  if (!spill.GetIsMC()) {
437 
438  header._POT_NoCut += beam.POTSincePreviousSavedSpill;
439  header._Spill_NoCut += beam.SpillsSincePreviousSavedSpill;
440 
441  if (beam.GoodSpill == 0) {
442  header._POT_BadBeam += beam.POTSincePreviousSavedSpill;
443  header._Spill_BadBeam += beam.SpillsSincePreviousSavedSpill;
444  return;
445  }
446 
447  if (!spill.DataQuality->GoodDaq) {
448  header._POT_BadND280 += beam.POTSincePreviousSavedSpill;
449  header._Spill_BadND280 += beam.SpillsSincePreviousSavedSpill;
450  return;
451  }
452 
453  if (beam.GoodSpill == 100) {
454  header._POT_0KA += beam.POTSincePreviousSavedSpill;
455  } else if (beam.GoodSpill == 1) {
456  header._POT_250KA += beam.POTSincePreviousSavedSpill;
457  } else if (beam.GoodSpill == 2) {
458  header._POT_200KA += beam.POTSincePreviousSavedSpill;
459  } else if (beam.GoodSpill == -1) {
460  header._POT_m250KA += beam.POTSincePreviousSavedSpill;
461  } else {
462  header._POT_OtherKA += beam.POTSincePreviousSavedSpill;
463  }
464 
465  header._POT_GoodBeamGoodND280 += beam.POTSincePreviousSavedSpill;
466  header._Spill_GoodBeamGoodND280 += beam.SpillsSincePreviousSavedSpill;
467  }
468  else{
469  // For MC there is no information about magnet current, Spill and DQ are always good
470  header._Spill_NoCut += beam.SpillsSincePreviousSavedSpill;
471  header._Spill_GoodBeamGoodND280 += beam.SpillsSincePreviousSavedSpill;
472 
473  header._POT_NoCut += beam.POTSincePreviousSavedSpill;
474  header._POT_GoodBeamGoodND280 += beam.POTSincePreviousSavedSpill;
475  }
476 
477 }
478 
479 //********************************************************************
480 Float_t anaUtils::ComputeInversePT(const AnaDetCrossingB &cross, bool entrance) {
481 //********************************************************************
482 
483  Float_t p = entrance ? GetEntranceMomentum(cross): GetExitMomentum(cross);
484 
485  if (p <= 0) return -999;
486 
487  Float_t ut = entrance ? sqrt(1-pow(cross.EntranceMomentum[0]/p,2)) :
488  sqrt(1-pow(cross.ExitMomentum[0]/p,2));
489 
490  Float_t pt = ut * p;
491 
492  if (pt!= 0)
493  return 1/pt;
494  else
495  return -999;
496 }
497 
498 //********************************************************************
499 Float_t anaUtils::ComputeInversePTFlip(const AnaTrackB& track){
500 //********************************************************************
501 
502  Float_t u_t = sqrt(1-pow(track.DirectionEnd[0],2));
503  Float_t pt = u_t*(track.MomentumFlip);
504  if (pt!=0)
505  return 1/pt;
506  else
507  return -999;
508 }
509 
510 //********************************************************************
511 Float_t anaUtils::ComputeMomentumFromInversePTFlip(const AnaParticleB &part, Float_t PTinv) {
512 //********************************************************************
513 
514  if (PTinv==0) return -999;
515  Float_t pt = 1/PTinv;
516  Float_t u_t = sqrt(1-pow(part.DirectionEnd[0],2));
517  if (u_t==0) return -999;
518 
519  Float_t p = pt/u_t;
520  return p;
521 }
522 
AnaBeamB
Representation of the beam quality and perhaps other beam information as needed.
Definition: BaseDataClasses.hxx:758
Header::_Spill_NoCut
Int_t _Spill_NoCut
Spill info.
Definition: Header.hxx:120
anaUtils::GetFGDMichelElectrons
int GetFGDMichelElectrons(const AnaEventB &event, const SubDetId::SubDetEnum det, AnaFgdTimeBinB **arr, bool prod5Cut=0)
Get all delayed time bins as Michel Electron candidates.
Definition: ND280AnalysisUtils.cxx:98
AnaParticleB::DirectionEnd
Float_t DirectionEnd[3]
The reconstructed end direction of the particle.
Definition: BaseDataClasses.hxx:289
anaUtils::GetExitMomentum
Float_t GetExitMomentum(const AnaDetCrossingB &cross)
Get the momentum value from AnaDetCrossing.
Definition: TruthUtils.cxx:459
AnaBeamB::SpillsSincePreviousSavedSpill
Int_t SpillsSincePreviousSavedSpill
Definition: BaseDataClasses.hxx:785
AnaSpillB::EventInfo
AnaEventInfoB * EventInfo
Run, sunrun, event, time stamp, etc.
Definition: BaseDataClasses.hxx:1022
Header
This class handles POT info, SoftwareVersion and IsMC.
Definition: Header.hxx:10
AnaEventInfoB::SubRun
Int_t SubRun
The subrun number.
Definition: BaseDataClasses.hxx:945
AnaTrackB::MomentumFlip
Float_t MomentumFlip
Momentum for the main PID hypothesis and reverse sense.
Definition: BaseDataClasses.hxx:638
anaUtils::GetEntranceMomentum
Float_t GetEntranceMomentum(const AnaDetCrossingB &cross)
Get the momentum value from AnaDetCrossing.
Definition: TruthUtils.cxx:451
anaUtils::ComputeInversePT
Float_t ComputeInversePT(const AnaDetCrossingB &cross, bool entrance=true)
Compute inverse PT given an AnaDetCrossing.
Definition: ND280AnalysisUtils.cxx:480
AnaSpillB::Beam
AnaBeamB * Beam
The beam quality flags for this spill.
Definition: BaseDataClasses.hxx:1028
anaUtils::GetRunPeriod
int GetRunPeriod(int run, int subrun=-1)
Returns the run period (sequentially: 0,1,2,3,4,5 ...)
Definition: ND280AnalysisUtils.cxx:17
anaUtils::ComputeMomentumFromInversePTFlip
Float_t ComputeMomentumFromInversePTFlip(const AnaParticleB &part, Float_t PTinv)
Compute the total momentum (flip) given the part and the inverse transverse momentum.
Definition: ND280AnalysisUtils.cxx:511
AnaDetCrossingB
Representation of a detector crossing info for a true particle (G4 trajectory).
Definition: BaseDataClasses.hxx:56
SubDetId::SubDetEnum
SubDetEnum
Enumeration of all detector systems and subdetectors.
Definition: SubDetId.hxx:25
BeamBunching::GetNBunches
UInt_t GetNBunches(const AnaEventB &event)
Number of bunches in the run period for the current run.
Definition: BeamBunching.cxx:61
AnaSpillB::GetIsMC
bool GetIsMC() const
Return whether this spill is from Monte Carlo or not.
Definition: BaseDataClasses.hxx:985
AnaTrackB
Representation of a global track.
Definition: BaseDataClasses.hxx:603
AnaDetCrossingB::ExitMomentum
Float_t ExitMomentum[3]
for each subdetector tell the exit momentum
Definition: BaseDataClasses.hxx:84
AnaBeamB::POTSincePreviousSavedSpill
Double_t POTSincePreviousSavedSpill
Definition: BaseDataClasses.hxx:781
AnaBeamB::GoodSpill
Int_t GoodSpill
Good spill flag, as defined in Beam Summary Data. 0 is bad.
Definition: BaseDataClasses.hxx:777
AnaEventInfoB::Event
Int_t Event
The ND280 event number.
Definition: BaseDataClasses.hxx:948
BeamBunching::GetBunchCentralTime
Float_t GetBunchCentralTime(const AnaEventB &event, Int_t ibunch)
Get the central time for bunch ibunch.
Definition: BeamBunching.cxx:77
anaUtils::GetSandMode
int GetSandMode(int run)
Returns the sans muon modes (sequentially: 0,1,2)
Definition: ND280AnalysisUtils.cxx:87
AnaParticleB
Representation of a reconstructed particle (track or shower).
Definition: BaseDataClasses.hxx:252
AnaSpillB::DataQuality
AnaDataQualityB * DataQuality
The ND280 data quality flags for this spill.
Definition: BaseDataClasses.hxx:1025
AnaDetCrossingB::EntranceMomentum
Float_t EntranceMomentum[3]
for each subdetector tell the entrance momentum
Definition: BaseDataClasses.hxx:81
AnaEventInfoB::Run
Int_t Run
The ND280 run number.
Definition: BaseDataClasses.hxx:942
Generated by   doxygen 1.8.13