1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | #include <math.h> |
13 | |
14 | #include "JEventProcessor_eta6g_skim.h" |
15 | using namespace jana; |
16 | |
17 | |
18 | #include "JANA/JApplication.h" |
19 | #include <TLorentzVector.h> |
20 | #include "TMath.h" |
21 | #include "JANA/JApplication.h" |
22 | #include "DANA/DApplication.h" |
23 | #include "FCAL/DFCALShower.h" |
24 | #include "BCAL/DBCALShower.h" |
25 | #include "CCAL/DCCALShower.h" |
26 | #include "FCAL/DFCALCluster.h" |
27 | #include "FCAL/DFCALHit.h" |
28 | #include "START_COUNTER/DSCHit.h" |
29 | #include "ANALYSIS/DAnalysisUtilities.h" |
30 | #include "PID/DVertex.h" |
31 | #include "PID/DEventRFBunch.h" |
32 | #include "HDDM/DEventWriterHDDM.h" |
33 | #include <HDGEOMETRY/DGeometry.h> |
34 | #include <TOF/DTOFPoint.h> |
35 | |
36 | #include "GlueX.h" |
37 | #include <vector> |
38 | #include <deque> |
39 | #include <string> |
40 | #include <iostream> |
41 | #include <algorithm> |
42 | #include <stdio.h> |
43 | #include <stdlib.h> |
44 | |
45 | extern "C"{ |
46 | void InitPlugin(JApplication *app){ |
47 | InitJANAPlugin(app); |
48 | app->AddProcessor(new JEventProcessor_eta6g_skim()); |
49 | } |
50 | } |
51 | |
52 | |
53 | |
54 | |
55 | |
56 | JEventProcessor_eta6g_skim::JEventProcessor_eta6g_skim() |
57 | { |
58 | |
59 | WRITE_EVIO = 1; |
60 | WRITE_HDDM = 0; |
61 | |
62 | gPARMS->SetDefaultParameter( "ETA6G_SKIM:WRITE_EVIO", WRITE_EVIO ); |
63 | gPARMS->SetDefaultParameter( "ETA6G_SKIM:WRITE_HDDM", WRITE_HDDM ); |
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 | JEventProcessor_eta6g_skim::~JEventProcessor_eta6g_skim() |
92 | { |
93 | |
94 | } |
95 | |
96 | |
97 | |
98 | |
99 | jerror_t JEventProcessor_eta6g_skim::init(void) |
100 | { |
101 | num_epics_events = 0; |
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 | combi6 = new Combination (6); |
128 | combi7 = new Combination (7); |
129 | |
130 | return NOERROR; |
131 | } |
132 | |
133 | |
134 | |
135 | |
136 | jerror_t JEventProcessor_eta6g_skim::brun(JEventLoop *eventLoop, int32_t runnumber) |
137 | { |
138 | DGeometry* dgeom = NULL__null; |
139 | DApplication* dapp = dynamic_cast< DApplication* >(eventLoop->GetJApplication()); |
140 | if (dapp) dgeom = dapp->GetDGeometry(runnumber); |
141 | if (dgeom) { |
142 | dgeom->GetTargetZ(m_targetZ); |
143 | } else { |
144 | cerr << "No geometry accessbile to ccal_timing monitoring plugin." << endl; |
145 | return RESOURCE_UNAVAILABLE; |
146 | } |
147 | jana::JCalibration *jcalib = japp->GetJCalibration(runnumber); |
148 | std::map<string, float> beam_spot; |
149 | jcalib->Get("PHOTON_BEAM/beam_spot", beam_spot); |
150 | m_beamSpotX = beam_spot.at("x"); |
151 | m_beamSpotY = beam_spot.at("y"); |
152 | |
153 | return NOERROR; |
154 | } |
155 | |
156 | |
157 | |
158 | |
159 | jerror_t JEventProcessor_eta6g_skim::evnt(JEventLoop *loop, uint64_t eventnumber) |
160 | { |
161 | |
162 | vector< const DFCALShower* > locFCALShowers; |
163 | vector< const DBCALShower* > locBCALShowers; |
164 | vector< const DCCALShower* > locCCALShowers; |
165 | vector< const DVertex* > kinfitVertex; |
166 | vector<const DTOFPoint*> tof_points; |
167 | vector< const DBeamPhoton* > locBeamPhotons; |
168 | vector< const DSCHit* > locSCHit; |
169 | vector<const DFCALHit *> locFCALHits; |
170 | |
171 | loop->Get(locFCALHits); |
172 | |
173 | loop->Get(locFCALShowers); |
174 | loop->Get(locBCALShowers); |
175 | loop->Get(locCCALShowers); |
176 | loop->Get(kinfitVertex); |
177 | loop->Get(tof_points); |
178 | loop->Get(locBeamPhotons); |
179 | loop->Get(locSCHit); |
180 | |
181 | vector<const DTrackCandidate*> locTrackCandidate; |
182 | loop->Get(locTrackCandidate); |
183 | |
184 | vector< const DTrackTimeBased* > locTrackTimeBased; |
185 | loop->Get(locTrackTimeBased); |
186 | |
187 | |
188 | |
189 | vector < const DFCALShower * > matchedShowers; |
190 | |
191 | const DEventWriterEVIO* locEventWriterEVIO = NULL__null; |
192 | loop->GetSingle(locEventWriterEVIO); |
193 | |
194 | |
195 | if(loop->GetJEvent().GetStatusBit(kSTATUS_BOR_EVENT)) { |
196 | |
197 | locEventWriterEVIO->Write_EVIOEvent( loop, "eta6g_skim" ); |
198 | return NOERROR; |
199 | } |
200 | |
201 | |
202 | if(loop->GetJEvent().GetStatusBit(kSTATUS_EPICS_EVENT) && (num_epics_events<5)) { |
203 | |
204 | locEventWriterEVIO->Write_EVIOEvent( loop, "eta6g_skim" ); |
205 | num_epics_events++; |
206 | return NOERROR; |
207 | } |
208 | |
209 | |
210 | vector< const JObject* > locObjectsToSave; |
211 | |
212 | bool Candidate = false; |
213 | |
214 | DVector3 vertex; |
215 | vertex.SetXYZ(m_beamSpotX, m_beamSpotY, m_targetZ); |
216 | for (unsigned int i = 0 ; i < tof_points.size(); i++) { |
217 | locObjectsToSave.push_back(static_cast<const JObject *>(tof_points[i])); |
218 | } |
219 | for (unsigned int i = 0 ; i < locBCALShowers.size(); i++) { |
220 | locObjectsToSave.push_back(static_cast<const JObject *>(locBCALShowers[i])); |
221 | } |
222 | for (unsigned int i = 0 ; i < locFCALShowers.size(); i++) { |
223 | locObjectsToSave.push_back(static_cast<const JObject *>(locFCALShowers[i])); |
224 | } |
225 | for (unsigned int i = 0 ; i < locCCALShowers.size(); i++) { |
226 | locObjectsToSave.push_back(static_cast<const JObject *>(locCCALShowers[i])); |
227 | } |
228 | for (unsigned int i = 0 ; i < locTrackCandidate.size(); i++) { |
229 | locObjectsToSave.push_back(static_cast<const JObject *>(locTrackCandidate[i])); |
230 | } |
231 | for (unsigned int i = 0 ; i < locTrackTimeBased.size(); i++) { |
232 | locObjectsToSave.push_back(static_cast<const JObject *>(locTrackTimeBased[i])); |
233 | } |
234 | for (unsigned int i = 0 ; i < locBeamPhotons.size(); i++) { |
235 | locObjectsToSave.push_back(static_cast<const JObject *>(locBeamPhotons[i])); |
236 | } |
237 | for (unsigned int i = 0 ; i < locSCHit.size(); i++) { |
238 | locObjectsToSave.push_back(static_cast<const JObject *>(locSCHit[i])); |
239 | } |
240 | |
241 | double kinfitVertexX = m_beamSpotX; |
242 | double kinfitVertexY = m_beamSpotY; |
243 | double kinfitVertexZ = m_targetZ; |
244 | for (unsigned int i = 0 ; i < kinfitVertex.size(); i++) { |
245 | if(i==0) |
246 | locObjectsToSave.push_back(static_cast<const JObject *>(kinfitVertex[0])); |
247 | } |
248 | |
249 | vector<const DEventRFBunch*> locEventRFBunches; |
250 | loop->Get(locEventRFBunches); |
251 | if(locEventRFBunches.size() > 0) { |
252 | locObjectsToSave.push_back(static_cast<const JObject *>(locEventRFBunches[0])); |
253 | } |
254 | |
255 | vector <TLorentzVector> PhotonList; PhotonList.clear(); |
256 | int photon_nb = locBCALShowers.size() + locFCALShowers.size() + locCCALShowers.size(); |
257 | for (unsigned int i = 0; i < locBCALShowers.size(); i ++) { |
258 | double e = locBCALShowers[i]->E; |
259 | double x = locBCALShowers[i]->x - kinfitVertexX; |
260 | double y = locBCALShowers[i]->y - kinfitVertexY; |
261 | double z = locBCALShowers[i]->z - kinfitVertexZ; |
262 | DVector3 vertex(x, y, z); |
263 | |
264 | |
265 | double p = e; |
266 | double px = p * sin(vertex.Theta()) * cos(vertex.Phi()); |
267 | double py = p * sin(vertex.Theta()) * sin(vertex.Phi()); |
268 | double pz = p * cos(vertex.Theta()); |
269 | TLorentzVector PhotonVec(px, py, pz, e); |
270 | if (e > 0.2) PhotonList.push_back(PhotonVec); |
271 | } |
272 | for (unsigned int i = 0; i < locFCALShowers.size(); i ++) { |
273 | double e = locFCALShowers[i]->getEnergy(); |
274 | double x = locFCALShowers[i]->getPosition().X() - kinfitVertexX; |
275 | double y = locFCALShowers[i]->getPosition().Y() - kinfitVertexY; |
276 | double z = locFCALShowers[i]->getPosition().Z() - kinfitVertexZ; |
277 | DVector3 vertex(x, y, z); |
278 | |
279 | |
280 | double p = e; |
281 | double px = p * sin(vertex.Theta()) * cos(vertex.Phi()); |
282 | double py = p * sin(vertex.Theta()) * sin(vertex.Phi()); |
283 | double pz = p * cos(vertex.Theta()); |
284 | TLorentzVector PhotonVec(px, py, pz, e); |
285 | if (e > 0.2) PhotonList.push_back(PhotonVec); |
286 | } |
287 | |
288 | |
289 | |
290 | |
291 | |
292 | |
293 | |
294 | |
295 | |
296 | |
297 | |
298 | |
299 | |
300 | |
301 | |
302 | |
303 | |
304 | Double_t bestChi2Eta = 1.0e30; |
305 | Double_t bestChi2EtaPrim = 1.0e30; |
306 | vector <TLorentzVector> PhotonEta6gList;PhotonEta6gList.clear(); |
307 | vector <TLorentzVector> PhotonEtaprim6gList;PhotonEtaprim6gList.clear(); |
308 | |
309 | |
310 | |
311 | |
312 | |
313 | |
314 | |
315 | |
316 | |
317 | |
318 | |
319 | |
320 | photon_nb = PhotonList.size(); |
321 | double FCAL_trg_Esum = 0; |
322 | |
323 | for (vector<const DFCALHit*>::const_iterator hit = locFCALHits.begin(); hit != locFCALHits.end(); hit++ ) { |
324 | if ((**hit).E > 0.150) |
325 | FCAL_trg_Esum += (**hit).E; |
326 | } |
327 | |
328 | |
329 | Candidate |= ( (6 <= photon_nb && photon_nb <= 15) ); |
| Value stored to 'Candidate' is never read |
330 | |
331 | if ( FCAL_trg_Esum > 3.5 && PhotonList.size() >= 6 ) { |
332 | |
333 | if( WRITE_EVIO ){ |
334 | |
335 | locEventWriterEVIO->Write_EVIOEvent( loop, "eta6g_skim"); |
336 | } |
337 | if( WRITE_HDDM ) { |
338 | vector<const DEventWriterHDDM*> locEventWriterHDDMVector; |
339 | loop->Get(locEventWriterHDDMVector); |
340 | locEventWriterHDDMVector[0]->Write_HDDMEvent(loop, ""); |
341 | } |
342 | |
343 | } |
344 | |
345 | |
346 | return NOERROR; |
347 | } |
348 | |
349 | |
350 | |
351 | |
352 | jerror_t JEventProcessor_eta6g_skim::erun(void) |
353 | { |
354 | |
355 | |
356 | |
357 | return NOERROR; |
358 | } |
359 | |
360 | |
361 | |
362 | |
363 | jerror_t JEventProcessor_eta6g_skim::fini(void) |
364 | { |
365 | |
366 | return NOERROR; |
367 | } |
368 | void JEventProcessor_eta6g_skim::Combined6g(vector<TLorentzVector>&EMList, |
369 | Double_t &bestChi0Eta, |
370 | Double_t &bestChi0EtaPrim, |
371 | vector<TLorentzVector>&PhotonEta6gList, |
372 | vector<TLorentzVector>&PhotonEtaprim6gList) |
373 | { |
374 | bestChi0EtaPrim = 1.0e30; |
375 | bestChi0Eta = 1.0e30; |
376 | if(EMList.size() == 6) { |
377 | for (int i_comb = 0; i_comb < (*combi6).getMaxCombi(); i_comb ++) { |
378 | combi6->getCombi(i_comb); |
379 | |
380 | double Esum = 0.0; |
381 | for (int i = 0; i < 6; i ++) { |
382 | Esum += EMList[combi6->combi[i]].E(); |
383 | } |
384 | |
385 | double Chi2_pi0Mass[3]; |
386 | double Chi2_etaMass[3]; |
387 | vector<TLorentzVector>GG;GG.clear(); |
388 | vector<TLorentzVector>Pi0Cor;Pi0Cor.clear(); |
389 | vector<TLorentzVector>EtaCor;EtaCor.clear(); |
390 | |
391 | for (int i = 0; i < 3; i ++) { |
392 | GG.push_back(EMList[combi6->combi[2*i]] + EMList[combi6->combi[2*i+1]] ); |
393 | Pi0Cor.push_back( pi0Mass / GG[i].M() * GG[i] ); |
394 | EtaCor.push_back( etaMass / GG[i].M() * GG[i] ); |
395 | Chi2_pi0Mass[i] = TMath::Power((GG[i].M() - pi0Mass) / 12.8e-3,2.0); |
396 | Chi2_etaMass[i] = TMath::Power((GG[i].M() - etaMass) / 31.1e-3,2.0); |
397 | } |
398 | |
399 | double Chi2_3pi0 = Chi2_pi0Mass[0] + Chi2_pi0Mass[1] + Chi2_pi0Mass[2]; |
400 | double Chi2_2pi0eta_0 = Chi2_pi0Mass[0] + Chi2_pi0Mass[1] + Chi2_etaMass[2]; |
401 | double Chi2_2pi0eta_1 = Chi2_pi0Mass[0] + Chi2_etaMass[1] + Chi2_pi0Mass[2]; |
402 | double Chi2_2pi0eta_2 = Chi2_etaMass[0] + Chi2_pi0Mass[1] + Chi2_pi0Mass[2]; |
403 | |
404 | if (Esum > 500.0e-3) { |
405 | TLorentzVector EtaVec = Pi0Cor[0] + Pi0Cor[1] + Pi0Cor[2]; |
406 | bool AnEta = false; |
407 | if (GG[0].M() > 110.0e-3 && |
408 | GG[1].M() > 110.0e-3 && |
409 | GG[2].M() > 110.0e-3 |
410 | && |
411 | GG[0].M() < 160.0e-3 && |
412 | GG[1].M() < 160.0e-3 && |
413 | GG[2].M() < 160.0e-3) |
414 | AnEta = true; |
415 | |
416 | bool AnEtaPrim = false; |
417 | bool AnEtaPrim1 = false; |
418 | bool AnEtaPrim2 = false; |
419 | bool AnEtaPrim3 = false; |
420 | if (GG[0].M() > 110.0e-3 && |
421 | GG[1].M() > 110.0e-3 && |
422 | GG[2].M() > 500.0e-3 |
423 | && |
424 | GG[0].M() < 160.0e-3 && |
425 | GG[1].M() < 160.0e-3 && |
426 | GG[2].M() < 600.0e-3) |
427 | AnEtaPrim1 = true; |
428 | if (GG[0].M() > 110.0e-3 && |
429 | GG[1].M() > 500.0e-3 && |
430 | GG[2].M() > 110.0e-3 |
431 | && |
432 | GG[0].M() < 160.0e-3 && |
433 | GG[1].M() < 600.0e-3 && |
434 | GG[2].M() < 160.0e-3) |
435 | AnEtaPrim2 = true; |
436 | if (GG[0].M() > 500.0e-3 && |
437 | GG[1].M() > 110.0e-3 && |
438 | GG[2].M() > 110.0e-3 |
439 | && |
440 | GG[0].M() < 600.0e-3 && |
441 | GG[1].M() < 160.0e-3 && |
442 | GG[2].M() < 160.0e-3) |
443 | AnEtaPrim3 = true; |
444 | if (AnEtaPrim1 || AnEtaPrim2 || AnEtaPrim3) |
445 | AnEtaPrim = true; |
446 | |
447 | if (AnEta && !AnEtaPrim && Esum > 500.0e-3) { |
448 | if (Chi2_3pi0<bestChi0Eta) { |
449 | bestChi0Eta = Chi2_3pi0; |
450 | PhotonEta6gList.clear(); |
451 | PhotonEta6gList.push_back(EMList[combi6->combi[0]]); |
452 | PhotonEta6gList.push_back(EMList[combi6->combi[1]]); |
453 | PhotonEta6gList.push_back(EMList[combi6->combi[2]]); |
454 | PhotonEta6gList.push_back(EMList[combi6->combi[3]]); |
455 | PhotonEta6gList.push_back(EMList[combi6->combi[4]]); |
456 | PhotonEta6gList.push_back(EMList[combi6->combi[5]]); |
457 | } |
458 | } |
459 | |
460 | if (!AnEta && AnEtaPrim && Esum > 900.0e-3) { |
461 | double Esum2gg1 = EMList[combi6->combi[4]].E() + EMList[combi6->combi[5]].E(); |
462 | if (AnEtaPrim1 && Esum2gg1 > 500.0e-3) { |
463 | if (Chi2_2pi0eta_0 < bestChi0EtaPrim) { |
464 | bestChi0EtaPrim = Chi2_2pi0eta_0; |
465 | PhotonEtaprim6gList.clear(); |
466 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[0]]); |
467 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[1]]); |
468 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[2]]); |
469 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[3]]); |
470 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[4]]); |
471 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[5]]); |
472 | } |
473 | } |
474 | if (AnEtaPrim2) { |
475 | double Esum2gg2 = EMList[combi6->combi[2]].E() + EMList[combi6->combi[3]].E(); |
476 | if (Chi2_2pi0eta_1 < bestChi0EtaPrim && Esum2gg2 > 500.0e-3) { |
477 | PhotonEtaprim6gList.clear(); |
478 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[0]]); |
479 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[1]]); |
480 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[4]]); |
481 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[5]]); |
482 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[2]]); |
483 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[3]]); |
484 | } |
485 | } |
486 | if (AnEtaPrim3) { |
487 | double Esum2gg3 = EMList[combi6->combi[0]].E() + EMList[combi6->combi[1]].E(); |
488 | if (Chi2_2pi0eta_2 < bestChi0EtaPrim && Esum2gg3 > 500.0e-3) { |
489 | PhotonEtaprim6gList.clear(); |
490 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[2]]); |
491 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[3]]); |
492 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[4]]); |
493 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[5]]); |
494 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[0]]); |
495 | PhotonEtaprim6gList.push_back(EMList[combi6->combi[1]]); |
496 | } |
497 | } |
498 | } |
499 | } |
500 | } |
501 | } |
502 | if(PhotonEta6gList.size()>0) |
503 | PhotonEtaprim6gList.clear(); |
504 | } |
505 | void JEventProcessor_eta6g_skim::Combined7g(vector<TLorentzVector>&EMList, |
506 | Double_t &bestChi0Eta, |
507 | Double_t &bestChi0EtaPrim, |
508 | vector<TLorentzVector>&PhotonEta6gList, |
509 | vector<TLorentzVector>&PhotonEtaprim6gList) |
510 | { |
511 | bestChi0EtaPrim = 1.0e30; |
512 | bestChi0Eta = 1.0e30; |
513 | if (EMList.size() == 7) { |
514 | for (int i_comb = 0; i_comb < (*combi7).getMaxCombi(); i_comb ++) { |
515 | combi7->getCombi(i_comb); |
516 | |
517 | double Esum = 0.0; |
518 | for (int i = 0; i < 6; i ++) { |
519 | Esum += EMList[combi7->combi[i]].E(); |
520 | } |
521 | |
522 | double Chi2_pi0Mass[3]; |
523 | double Chi2_etaMass[3]; |
524 | vector<TLorentzVector>GG;GG.clear(); |
525 | vector<TLorentzVector>Pi0Cor;Pi0Cor.clear(); |
526 | vector<TLorentzVector>EtaCor;EtaCor.clear(); |
527 | |
528 | for (int i = 0; i < 3; i ++) { |
529 | GG.push_back( EMList[combi7->combi[2*i]] + EMList[combi7->combi[2*i+1]] ); |
530 | Pi0Cor.push_back( pi0Mass / GG[i].M() * GG[i] ); |
531 | EtaCor.push_back( etaMass / GG[i].M() * GG[i] ); |
532 | Chi2_pi0Mass[i] = TMath::Power((GG[i].M() - pi0Mass) / 12.8e-3,2.0); |
533 | Chi2_etaMass[i] = TMath::Power((GG[i].M() - etaMass) / 31.1e-3,2.0); |
534 | } |
535 | |
536 | if (Esum > 500.0e-3) { |
537 | TLorentzVector EtaVec = Pi0Cor[0] + Pi0Cor[1] + Pi0Cor[2]; |
538 | TLorentzVector EtaprimVec1 = Pi0Cor[0] + Pi0Cor[1] + EtaCor[2]; |
539 | TLorentzVector EtaprimVec2 = Pi0Cor[0] + EtaCor[1] + Pi0Cor[2]; |
540 | TLorentzVector EtaprimVec3 = EtaCor[0] + Pi0Cor[1] + Pi0Cor[2]; |
541 | TLorentzVector Candidate = (TLorentzVector) EMList[combi7->combi[6]]; |
542 | Double_t CopAngleEta = (fabs(EtaVec.Phi() - Candidate.Phi())) * TMath::RadToDeg(); |
543 | Double_t Chi2AngleEta = TMath::Power((180.0 - CopAngleEta) / 30.0,2.0); |
544 | |
545 | Double_t CopAngleEtaprim1 = (fabs(EtaprimVec1.Phi() - Candidate.Phi())) * TMath::RadToDeg(); |
546 | Double_t Chi2AngleEtaprim1 = TMath::Power((180.0 - CopAngleEtaprim1) / 30.0,2.0); |
547 | |
548 | Double_t CopAngleEtaprim2 = (fabs(EtaprimVec2.Phi() - Candidate.Phi())) * TMath::RadToDeg(); |
549 | Double_t Chi2AngleEtaprim2 = TMath::Power((180.0 - CopAngleEtaprim2) / 30.0,2.0); |
550 | |
551 | Double_t CopAngleEtaprim3 = (fabs(EtaprimVec3.Phi() - Candidate.Phi())) * TMath::RadToDeg(); |
552 | Double_t Chi2AngleEtaprim3 = TMath::Power((180.0 - CopAngleEtaprim3) / 30.0,2.0); |
553 | |
554 | double Chi2_3pi0 = Chi2_pi0Mass[0] + Chi2_pi0Mass[1] + Chi2_pi0Mass[2] + Chi2AngleEta; |
555 | double Chi2_2pi0eta_0 = Chi2_pi0Mass[0] + Chi2_pi0Mass[1] + Chi2_etaMass[2] + Chi2AngleEtaprim1; |
556 | double Chi2_2pi0eta_1 = Chi2_pi0Mass[0] + Chi2_etaMass[1] + Chi2_pi0Mass[2] + Chi2AngleEtaprim2; |
557 | double Chi2_2pi0eta_2 = Chi2_etaMass[0] + Chi2_pi0Mass[1] + Chi2_pi0Mass[2] + Chi2AngleEtaprim3; |
558 | |
559 | |
560 | |
561 | |
562 | |
563 | |
564 | bool AnEta = false; |
565 | if (GG[0].M() > 110.0e-3 && |
566 | GG[1].M() > 110.0e-3 && |
567 | GG[2].M() > 110.0e-3 |
568 | && |
569 | GG[0].M() < 160.0e-3 && |
570 | GG[1].M() < 160.0e-3 && |
571 | GG[2].M() < 160.0e-3) |
572 | AnEta = true; |
573 | |
574 | bool AnEtaPrim = false; |
575 | bool AnEtaPrim1 = false; |
576 | bool AnEtaPrim2 = false; |
577 | bool AnEtaPrim3 = false; |
578 | if (GG[0].M() > 110.0e-3 && |
579 | GG[1].M() > 110.0e-3 && |
580 | GG[2].M() > 500.0e-3 |
581 | && |
582 | GG[0].M() < 160.0e-3 && |
583 | GG[1].M() < 160.0e-3 && |
584 | GG[2].M() < 600.0e-3) |
585 | AnEtaPrim1 = true; |
586 | if (GG[0].M() > 110.0e-3 && |
587 | GG[1].M() > 500.0e-3 && |
588 | GG[2].M() > 110.0e-3 |
589 | && |
590 | GG[0].M() < 160.0e-3 && |
591 | GG[1].M() < 600.0e-3 && |
592 | GG[2].M() < 160.0e-3) |
593 | AnEtaPrim2 = true; |
594 | if (GG[0].M() > 500.0e-3 && |
595 | GG[1].M() > 110.0e-3 && |
596 | GG[2].M() > 110.0e-3 |
597 | && |
598 | GG[0].M() < 600.0e-3 && |
599 | GG[1].M() < 160.0e-3 && |
600 | GG[2].M() < 160.0e-3) |
601 | AnEtaPrim3 = true; |
602 | if(AnEtaPrim1 || AnEtaPrim2 || AnEtaPrim3) |
603 | AnEtaPrim = true; |
604 | |
605 | if (AnEta && !AnEtaPrim && Esum > 500.0e-3) { |
606 | if (Chi2_3pi0 < bestChi0Eta) { |
607 | bestChi0Eta = Chi2_3pi0; |
608 | PhotonEta6gList.clear(); |
609 | PhotonEta6gList.push_back(EMList[combi7->combi[0]]); |
610 | PhotonEta6gList.push_back(EMList[combi7->combi[1]]); |
611 | PhotonEta6gList.push_back(EMList[combi7->combi[2]]); |
612 | PhotonEta6gList.push_back(EMList[combi7->combi[3]]); |
613 | PhotonEta6gList.push_back(EMList[combi7->combi[4]]); |
614 | PhotonEta6gList.push_back(EMList[combi7->combi[5]]); |
615 | PhotonEta6gList.push_back(EMList[combi7->combi[6]]); |
616 | } |
617 | } |
618 | |
619 | if (!AnEta && AnEtaPrim && Esum > 900.0e-3) { |
620 | double Esum2gg1 = EMList[combi7->combi[4]].E() + EMList[combi7->combi[5]].E(); |
621 | if (AnEtaPrim1 && Esum2gg1 > 500.0e-3) { |
622 | if(Chi2_2pi0eta_0 < bestChi0EtaPrim) { |
623 | bestChi0EtaPrim = Chi2_2pi0eta_0; |
624 | PhotonEtaprim6gList.clear(); |
625 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[0]]); |
626 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[1]]); |
627 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[2]]); |
628 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[3]]); |
629 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[4]]); |
630 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[5]]); |
631 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[6]]); |
632 | } |
633 | } |
634 | double Esum2gg2 = EMList[combi7->combi[2]].E() + EMList[combi7->combi[3]].E(); |
635 | if (AnEtaPrim2 && Esum2gg2 > 500.0e-3) { |
636 | if(Chi2_2pi0eta_1 < bestChi0EtaPrim) { |
637 | PhotonEtaprim6gList.clear(); |
638 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[0]]); |
639 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[1]]); |
640 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[4]]); |
641 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[5]]); |
642 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[2]]); |
643 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[3]]); |
644 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[6]]); |
645 | } |
646 | } |
647 | double Esum2gg3 = EMList[combi7->combi[0]].E() + EMList[combi7->combi[1]].E(); |
648 | if (AnEtaPrim3 && Esum2gg3 > 500.0e-3) { |
649 | if(Chi2_2pi0eta_2 < bestChi0EtaPrim) { |
650 | PhotonEtaprim6gList.clear(); |
651 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[2]]); |
652 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[3]]); |
653 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[4]]); |
654 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[5]]); |
655 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[0]]); |
656 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[1]]); |
657 | PhotonEtaprim6gList.push_back(EMList[combi7->combi[6]]); |
658 | } |
659 | } |
660 | } |
661 | } |
662 | } |
663 | } |
664 | if(PhotonEta6gList.size()>0) |
665 | PhotonEtaprim6gList.clear(); |
666 | } |