BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
h3lAna.C
Go to the documentation of this file.
1
2#include "TFile.h"
3#include "TH1F.h"
4#include "TTreeReader.h"
5#include "TTreeReaderValue.h"
6
7#include <TLorentzVector.h>
8#include <TSystem.h>
9// #include "BmnHypNuclPair.h"
10
11using namespace std;
12using namespace TMath;
13
14Float_t invMassMin = 2.95;
15Float_t invMassMax = 3.05;
16
17const Int_t nPath = 1;
18const Int_t nDca2 = 1;
19const Int_t nDca12 = 1;
20const Int_t nDca0 = 1;
21
22// EXP max by S
23// vector<Float_t> pathMin = {1.0};
24// vector<Float_t> pathMax = {35};
25// vector<Float_t> dca2Min = {0.0};
26// vector<Float_t> dca2Max = {100};
27// vector<Float_t> dca12Min = {0.0};
28// vector<Float_t> dca12Max = {9};
29// vector<Float_t> dca0Min = {0.08};
30// vector<Float_t> dca0Max = {0.65};
31
32// about 58% of max S and 19% of max B (SIGN = 8.49, dSIGN = 1.37, S = 1262, dS = 207, B = 20854, dB = 144)
33// vector<Float_t> pathMin = {1.4};
34// vector<Float_t> pathMax = {35};
35// vector<Float_t> dca2Min = {3.6};
36// vector<Float_t> dca2Max = {100};
37// vector<Float_t> dca12Min = {0.0};
38// vector<Float_t> dca12Max = {6.0};
39// vector<Float_t> dca0Min = {0.08};
40// vector<Float_t> dca0Max = {0.65};
41
42// MC max S
43// vector<Float_t> pathMin = {1.2};
44// vector<Float_t> pathMax = {35};
45// vector<Float_t> dca2Min = {0.4};
46// vector<Float_t> dca2Max = {100};
47// vector<Float_t> dca12Min = {0.0};
48// vector<Float_t> dca12Max = {18};
49// vector<Float_t> dca0Min = {0.08};
50// vector<Float_t> dca0Max = {1.05};
51
52// MC by EXP (projection through param scan) (SIGN = 7.8, dSIGN = 1.1, S = 136, dS = 22, B = 164, dB = 13)
53// vector<Float_t> pathMin = {1.1};
54// vector<Float_t> pathMax = {35};
55// vector<Float_t> dca2Min = {3.6};
56// vector<Float_t> dca2Max = {100};
57// vector<Float_t> dca12Min = {0.0};
58// vector<Float_t> dca12Max = {6.0};
59// vector<Float_t> dca0Min = {0.08};
60// vector<Float_t> dca0Max = {0.65};
61
62vector<Float_t> pathMin = {0.3};
63vector<Float_t> pathMax = {35};
64vector<Float_t> dca2Min = {5.7};
65vector<Float_t> dca2Max = {14};
66vector<Float_t> dca12Min = {0.0};
67vector<Float_t> dca12Max = {5.0};
68vector<Float_t> dca0Min = {0.11};
69vector<Float_t> dca0Max = {2.5};
70
71const Int_t nMom1 = 1;
72vector<Float_t> mom1Min = {4.5};
73vector<Float_t> mom1Max = {15.5};
74const Int_t nMom2 = 1;
75vector<Float_t> mom2Min = {0.3};
76vector<Float_t> mom2Max = {1.1};
77
78const Int_t nFsd1 = 1;
79vector<Int_t> fsd1Min = {0};
80vector<Int_t> fsd1Max = {4};
81const Int_t nFsd2 = 1;
82vector<Int_t> fsd2Min = {0};
83vector<Int_t> fsd2Max = {4};
84const Int_t nGem1 = 1;
85vector<Int_t> gem1Min = {5};
86vector<Int_t> gem1Max = {7};
87const Int_t nGem2 = 1;
88vector<Int_t> gem2Min = {1};
89vector<Int_t> gem2Max = {7};
90
91void h3lAna(TString inputFileName = "/home/merz/data/h3l_sim/h3lPairs_XeCsI_3.8AGeV_mb_4mln.root",
92 TString outputFileName = "/home/merz/data/h3l_sim/h3lPlots_XeCsI_3.8AGeV_mb_4mln.root")
93{
94
95 if (gSystem->AccessPathName(inputFileName))
96 return;
97
98 TH1F* hInvMass[nPath][nDca2][nDca12][nDca0][nMom1][nMom2][nFsd1][nFsd2][nGem1][nGem2];
99 for (Int_t iPath = 0; iPath < nPath; iPath++) {
100 for (Int_t iDca2 = 0; iDca2 < nDca2; iDca2++) {
101 for (Int_t iDca12 = 0; iDca12 < nDca12; iDca12++) {
102 for (Int_t iDca0 = 0; iDca0 < nDca0; iDca0++) {
103 for (Int_t iMom1 = 0; iMom1 < nMom1; iMom1++) {
104 for (Int_t iMom2 = 0; iMom2 < nMom2; iMom2++) {
105 for (Int_t iFsd1 = 0; iFsd1 < nFsd1; iFsd1++) {
106 for (Int_t iFsd2 = 0; iFsd2 < nFsd2; iFsd2++) {
107 for (Int_t iGem1 = 0; iGem1 < nGem1; iGem1++) {
108 for (Int_t iGem2 = 0; iGem2 < nGem2; iGem2++) {
109 hInvMass[iPath][iDca2][iDca12][iDca0][iMom1][iMom2][iFsd1][iFsd2][iGem1]
110 [iGem2] = new TH1F(Form("hInvMass_%d_%d_%d_%d_%d_%d_%d_%d_%d_%d",
111 iPath, iDca2, iDca12, iDca0, iMom1, iMom2,
112 iFsd1, iFsd2, iGem1, iGem2),
113 "; M; Count", 100, invMassMin, invMassMax);
114 }
115 }
116 }
117 }
118 }
119 }
120 }
121 }
122 }
123 }
124
125 cout << "File: " << inputFileName << endl;
126 TFile* myFile = TFile::Open(inputFileName);
127
128 TTreeReader myReader("bmndata", myFile);
129 TTreeReaderArray<BmnHypNuclPair> pairs(myReader, "pairs");
130
131 int evId = 0;
132
133 while (myReader.Next()) {
134 if (evId % 50000 == 0)
135 printf("evId = %d\n", evId);
136 evId++;
137 for (BmnHypNuclPair pair : pairs) {
138
139 // TLorentzVector lPos, lNeg;
140 // Float_t px1 = pair.GetPx1();
141 // Float_t py1 = pair.GetPy1();
142 // Float_t pz1 = pair.GetPz1();
143 // Float_t tx1 = px1 / pz1;
144 // Float_t ty1 = py1 / pz1;
145 // Float_t mom1_new = mom1; // - 50.0 * corrFit->Eval(mom1);
146 // //cout << mom1 << "\t" << corrFit->Eval(mom1) << "\t" << mom1_new << endl;
147 // Float_t pz1_new = mom1_new / TMath::Sqrt(tx1 * tx1 + ty1 * ty1 + 1);
148 // Float_t px1_new = pz1_new * tx1;
149 // Float_t py1_new = pz1_new * ty1;
150 // Float_t px2 = pair.GetPx2();
151 // Float_t py2 = pair.GetPy2();
152 // Float_t pz2 = pair.GetPz2();
153
154 // // if ((px1 * px2) > 0 || (py1 * py2) > 0)
155 // // continue;
156
157 // //cout << px1 << "\t" << py1 << "\t" << pz1 << endl;
158 // //cout << px2 << "\t" << py2 << "\t" << pz2 << endl;
159 // //cout << px1_new << "\t" << py1_new << "\t" << pz1_new << endl;
160 // lPos.SetXYZM(px1, py1, pz1, 2.8094135);
161 // lNeg.SetXYZM(px2, py2, pz2, 0.1395704);
162 // TLorentzVector lHyp = lPos + lNeg;
163 // Float_t pt = lHyp.Pt();
164
165 // if (pt < 0.6)
166 // continue;
167
168 // Float_t mass_new = TLorentzVector((lPos + lNeg)).Mag();
169 Float_t mass = pair.GetInvMass();
170 // Float_t mass_new = mass;
171 // cout << mass << "\t" << mass_new - mass << "\t" << mass_new << endl;
172 if (mass > invMassMax)
173 continue;
174 if (mass < invMassMin)
175 continue;
176 Float_t mom1 = pair.GetMom1();
177 Float_t mom2 = pair.GetMom2();
178 // Float_t massTof = pair.GetM2_700_1();
179 // if (massTof < m2CutLow->Eval(mom1_new) || massTof > m2CutUp->Eval(mom1_new))
180 // continue;
181 Float_t path = pair.GetPath();
182 Float_t dca0 = pair.GetDCA0();
183 Float_t dca12 = pair.GetDCA12() / pair.GetdDCA12();
184 Float_t dca2 = pair.GetDCA2() / pair.GetdDCA2();
185 Float_t dca1 = pair.GetDCA1() / pair.GetdDCA1();
186
187 Int_t fsd1 = pair.GetNHitsPart1("SILICON");
188 Int_t fsd2 = pair.GetNHitsPart2("SILICON");
189 Int_t gem1 = pair.GetNHitsPart1("GEM");
190 Int_t gem2 = pair.GetNHitsPart2("GEM");
191
192 for (Int_t iPath = 0; iPath < nPath; iPath++) {
193 if (path < pathMin[iPath] || path > pathMax[iPath])
194 continue;
195 for (Int_t iDca2 = 0; iDca2 < nDca2; iDca2++) {
196 if (dca2 < dca2Min[iDca2] || dca2 > dca2Max[iDca2])
197 continue;
198 for (Int_t iDca12 = 0; iDca12 < nDca12; iDca12++) {
199 if (dca12 < dca12Min[iDca12] || dca12 > dca12Max[iDca12])
200 continue;
201 for (Int_t iDca0 = 0; iDca0 < nDca0; iDca0++) {
202 if (dca0 < dca0Min[iDca0] || dca0 > dca0Max[iDca0])
203 continue;
204 for (Int_t iMom1 = 0; iMom1 < nMom1; iMom1++) {
205 if (mom1 < mom1Min[iMom1] || mom1 > mom1Max[iMom1])
206 continue;
207 for (Int_t iMom2 = 0; iMom2 < nMom2; iMom2++) {
208 if (mom2 < mom2Min[iMom2] || mom2 > mom2Max[iMom2])
209 continue;
210 for (Int_t iFsd1 = 0; iFsd1 < nFsd1; iFsd1++) {
211 if (fsd1 < fsd1Min[iFsd1] || fsd1 > fsd1Max[iFsd1])
212 continue;
213 for (Int_t iFsd2 = 0; iFsd2 < nFsd2; iFsd2++) {
214 if (fsd2 < fsd2Min[iFsd2] || fsd2 > fsd2Max[iFsd2])
215 continue;
216 for (Int_t iGem1 = 0; iGem1 < nGem1; iGem1++) {
217 if (gem1 < gem1Min[iGem1] || gem1 > gem1Max[iGem1])
218 continue;
219 for (Int_t iGem2 = 0; iGem2 < nGem2; iGem2++) {
220 if (gem2 < gem2Min[iGem2] || gem2 > gem2Max[iGem2])
221 continue;
222 hInvMass[iPath][iDca2][iDca12][iDca0][iMom1][iMom2][iFsd1][iFsd2]
223 [iGem1][iGem2]
224 ->Fill(mass);
225 }
226 }
227 }
228 }
229 }
230 }
231 }
232 }
233 }
234 }
235 }
236 }
237
238 myFile->Close();
239
240 TFile outFile(outputFileName, "RECREATE");
241 for (Int_t iPath = 0; iPath < nPath; iPath++) {
242 for (Int_t iDca2 = 0; iDca2 < nDca2; iDca2++) {
243 for (Int_t iDca12 = 0; iDca12 < nDca12; iDca12++) {
244 for (Int_t iDca0 = 0; iDca0 < nDca0; iDca0++) {
245 for (Int_t iMom1 = 0; iMom1 < nMom1; iMom1++) {
246 for (Int_t iMom2 = 0; iMom2 < nMom2; iMom2++) {
247 for (Int_t iFsd1 = 0; iFsd1 < nFsd1; iFsd1++) {
248 for (Int_t iFsd2 = 0; iFsd2 < nFsd2; iFsd2++) {
249 for (Int_t iGem1 = 0; iGem1 < nGem1; iGem1++) {
250 for (Int_t iGem2 = 0; iGem2 < nGem2; iGem2++) {
251 hInvMass[iPath][iDca2][iDca12][iDca0][iMom1][iMom2][iFsd1][iFsd2][iGem1]
252 [iGem2]
253 ->Write();
254 }
255 }
256 }
257 }
258 }
259 }
260 }
261 }
262 }
263 }
264 outFile.Close();
265}
Int_t
int Int_t
Float_t
float Float_t
Form
char * Form(const char *fmt,...)
gSystem
R__EXTERN TSystem * gSystem
TFile::Open
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
TFile::Close
void Close(Option_t *option="") override
TH1::Fill
virtual Int_t Fill(const char *name, Double_t w)
TObject::Write
virtual Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0)
nDca12
const Int_t nDca12
Definition h3lAna.C:19
nPath
const Int_t nPath
Definition h3lAna.C:17
nGem1
const Int_t nGem1
Definition h3lAna.C:84
nFsd2
const Int_t nFsd2
Definition h3lAna.C:81
nDca2
const Int_t nDca2
Definition h3lAna.C:18
nMom1
const Int_t nMom1
Definition h3lAna.C:71
nMom2
const Int_t nMom2
Definition h3lAna.C:74
nGem2
const Int_t nGem2
Definition h3lAna.C:87
nDca0
const Int_t nDca0
Definition h3lAna.C:20
nFsd1
const Int_t nFsd1
Definition h3lAna.C:78
invMassMax
Float_t invMassMax
Definition h3lAna.C:15
invMassMin
Float_t invMassMin
Definition h3lAna.C:14
h3lAna
void h3lAna(TString inputFileName="/home/merz/data/h3l_sim/h3lPairs_XeCsI_3.8AGeV_mb_4mln.root", TString outputFileName="/home/merz/data/h3l_sim/h3lPlots_XeCsI_3.8AGeV_mb_4mln.root")
Definition h3lAna.C:91
std
Definition json.hpp:4446