Bug Summary

File:plugins/Utilities/eta6g_primexd_skim/JEventProcessor_eta6g_skim.cc
Location:line 329, column 3
Description:Value stored to 'Candidate' is never read

Annotated Source Code

1/**************************************************************************
2* HallD software *
3* Copyright(C) 2020 GlueX and PrimEX-D Collaborations *
4* *
5* Author: The GlueX and PrimEX-D Collaborations *
6* Contributors: Igal Jaegle *
7* *
8* *
9* This software is provided "as is" without any warranty. *
10**************************************************************************/
11
12#include <math.h>
13
14#include "JEventProcessor_eta6g_skim.h"
15using namespace jana;
16
17// Routine used to create our JEventProcessor
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
45extern "C"{
46 void InitPlugin(JApplication *app){
47 InitJANAPlugin(app);
48 app->AddProcessor(new JEventProcessor_eta6g_skim());
49 }
50} // "C"
51
52
53//------------------
54// JEventProcessor_eta6g_skim (Constructor)
55//------------------
56JEventProcessor_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 MIN_MASS = 0.03; // GeV
67 MAX_MASS = 0.30; // GeV
68 MIN_E = 1.0; // GeV (photon energy cut)
69 MIN_R = 20; // cm (cluster distance to beam line)
70 MAX_DT = 10; // ns (cluster time diff. cut)
71 MAX_ETOT = 12; // GeV (max total FCAL energy)
72 MIN_BLOCKS = 2; // minumum blocks per cluster
73
74 WRITE_ROOT = 0;
75 WRITE_EVIO = 1;
76
77 gPARMS->SetDefaultParameter( "ETA6G_SKIM:MIN_MASS", MIN_MASS );
78 gPARMS->SetDefaultParameter( "ETA6G_SKIM:MAX_MASS", MAX_MASS );
79 gPARMS->SetDefaultParameter( "ETA6G_SKIM:MIN_E", MIN_E );
80 gPARMS->SetDefaultParameter( "ETA6G_SKIM:MIN_R", MIN_R );
81 gPARMS->SetDefaultParameter( "ETA6G_SKIM:MAX_DT", MAX_DT );
82 gPARMS->SetDefaultParameter( "ETA6G_SKIM:MAX_ETOT", MAX_ETOT );
83 gPARMS->SetDefaultParameter( "ETA6G_SKIM:MIN_BLOCKS", MIN_BLOCKS );
84 gPARMS->SetDefaultParameter( "ETA6G_SKIM:WRITE_ROOT", WRITE_ROOT );
85 */
86}
87
88//------------------
89// ~JEventProcessor_eta6g_skim (Destructor)
90//------------------
91JEventProcessor_eta6g_skim::~JEventProcessor_eta6g_skim()
92{
93
94}
95
96//------------------
97// init
98//------------------
99jerror_t JEventProcessor_eta6g_skim::init(void)
100{
101 num_epics_events = 0;
102/*
103 if( ! ( WRITE_ROOT || WRITE_EVIO ) ){
104
105 cerr << "No output mechanism has been specified." << endl;
106 return UNRECOVERABLE_ERROR;
107 }
108
109 if( WRITE_ROOT ){
110
111 japp->RootWriteLock();
112
113 m_tree = new TTree( "cluster", "Cluster Tree for Pi0 Calibration" );
114 m_tree->Branch( "nClus", &m_nClus, "nClus/I" );
115 m_tree->Branch( "hit0", m_hit0, "hit0[nClus]/I" );
116 m_tree->Branch( "px", m_px, "px[nClus]/F" );
117 m_tree->Branch( "py", m_py, "py[nClus]/F" );
118 m_tree->Branch( "pz", m_pz, "pz[nClus]/F" );
119
120 m_tree->Branch( "nHit", &m_nHit, "nHit/I" );
121 m_tree->Branch( "chan", m_chan, "chan[nHit]/I" );
122 m_tree->Branch( "e", m_e, "e[nHit]/F" );
123
124 japp->RootUnLock();
125 }
126*/
127 combi6 = new Combination (6);
128 combi7 = new Combination (7);
129
130 return NOERROR;
131}
132
133//------------------
134// brun
135//------------------
136jerror_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// evnt
158//------------------
159jerror_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 // always write out BOR events
195 if(loop->GetJEvent().GetStatusBit(kSTATUS_BOR_EVENT)) {
196 //jout << "Found BOR!" << endl;
197 locEventWriterEVIO->Write_EVIOEvent( loop, "eta6g_skim" );
198 return NOERROR;
199 }
200
201 // write out the first few EPICS events to save run number & other meta info
202 if(loop->GetJEvent().GetStatusBit(kSTATUS_EPICS_EVENT) && (num_epics_events<5)) {
203 //jout << "Found EPICS!" << endl;
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 //Use kinfit when available
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 //double r = vertex.Mag();
264 //double t = locCCALShowers[j]->time - (r / TMath::C());
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 //double r = vertex1.Mag();
279 //double t = locFCALShowers[i]->getTime() - (r / TMath::C());
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 for (unsigned int i = 0; i < locCCALShowers.size(); i ++) {
289 double e = locCCALShowers[i]->E;
290 double x = locCCALShowers[i]->x - kinfitVertexX;
291 double y = locCCALShowers[i]->y - kinfitVertexY;
292 double z = locCCALShowers[i]->z - kinfitVertexZ;
293 DVector3 vertex(x, y, z);
294 //double r = vertex.Mag();
295 //double t = locCCALShowers[j]->time - (r / TMath::C());
296 double p = e;
297 double px = p * sin(vertex.Theta()) * cos(vertex.Phi());
298 double py = p * sin(vertex.Theta()) * sin(vertex.Phi());
299 double pz = p * cos(vertex.Theta());
300 TLorentzVector PhotonVec(px, py, pz, e);
301 PhotonList.push_back(PhotonVec);
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 Combined6g(PhotonList,
310 bestChi2Eta,
311 bestChi2EtaPrim,
312 PhotonEta6gList,
313 PhotonEtaprim6gList);
314 Combined7g(PhotonList,
315 bestChi2Eta,
316 bestChi2EtaPrim,
317 PhotonEta6gList,
318 PhotonEtaprim6gList);
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 //cout <<"eta6g_skim"<<endl;
333 if( WRITE_EVIO ){
334 //locEventWriterEVIO->Write_EVIOEvent( loop, "eta6g_skim", locObjectsToSave );
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// erun
351//------------------
352jerror_t JEventProcessor_eta6g_skim::erun(void)
353{
354 // This is called whenever the run number changes, before it is
355 // changed to give you a chance to clean up before processing
356 // events from the next run number.
357 return NOERROR;
358}
359
360//------------------
361// Fin
362//------------------
363jerror_t JEventProcessor_eta6g_skim::fini(void)
364{
365 // Called before program exit after event processing is finished.
366 return NOERROR;
367}
368void 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}
505void 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 double Chi2_3pi0 = Chi2_pi0Mass[0] + Chi2_pi0Mass[1] + Chi2_pi0Mass[2];
560 double Chi2_2pi0eta_0 = Chi2_pi0Mass[0] + Chi2_pi0Mass[1] + Chi2_etaMass[2];
561 double Chi2_2pi0eta_1 = Chi2_pi0Mass[0] + Chi2_etaMass[1] + Chi2_pi0Mass[2];
562 double Chi2_2pi0eta_2 = Chi2_etaMass[0] + Chi2_pi0Mass[1] + Chi2_pi0Mass[2];
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}