Bug Summary

File:plugins/Utilities/pi0fcalskim/JEventProcessor_pi0fcalskim.cc
Location:line 200, column 4
Description:Value stored to 'p' is never read

Annotated Source Code

1// $Id$
2//
3// File: JEventProcessor_pi0fcalskim.cc
4// Created: Mon Dec 1 14:57:11 EST 2014
5// Creator: shepherd (on Linux ifarm1101 2.6.32-220.7.1.el6.x86_64 x86_64)
6//
7
8#include <math.h>
9
10#include "JEventProcessor_pi0fcalskim.h"
11using namespace jana;
12
13// Routine used to create our JEventProcessor
14#include "JANA/JApplication.h"
15#include <TLorentzVector.h>
16#include "TMath.h"
17#include "JANA/JApplication.h"
18#include "DANA/DApplication.h"
19#include "FCAL/DFCALShower.h"
20#include "FCAL/DFCALCluster.h"
21#include "FCAL/DFCALHit.h"
22#include "ANALYSIS/DAnalysisUtilities.h"
23#include "PID/DVertex.h"
24#include "GlueX.h"
25#include <vector>
26#include <deque>
27#include <string>
28#include <iostream>
29#include <algorithm>
30#include <stdio.h>
31#include <stdlib.h>
32
33extern "C"{
34 void InitPlugin(JApplication *app){
35 InitJANAPlugin(app);
36 app->AddProcessor(new JEventProcessor_pi0fcalskim());
37 }
38} // "C"
39
40
41//------------------
42// JEventProcessor_pi0fcalskim (Constructor)
43//------------------
44JEventProcessor_pi0fcalskim::JEventProcessor_pi0fcalskim()
45{
46
47 WRITE_EVIO = 1;
48/*
49 MIN_MASS = 0.03; // GeV
50 MAX_MASS = 0.30; // GeV
51 MIN_E = 1.0; // GeV (photon energy cut)
52 MIN_R = 20; // cm (cluster distance to beam line)
53 MAX_DT = 10; // ns (cluster time diff. cut)
54 MAX_ETOT = 12; // GeV (max total FCAL energy)
55 MIN_BLOCKS = 2; // minumum blocks per cluster
56
57 WRITE_ROOT = 0;
58 WRITE_EVIO = 1;
59
60 gPARMS->SetDefaultParameter( "PI0FCALSKIM:MIN_MASS", MIN_MASS );
61 gPARMS->SetDefaultParameter( "PI0FCALSKIM:MAX_MASS", MAX_MASS );
62 gPARMS->SetDefaultParameter( "PI0FCALSKIM:MIN_E", MIN_E );
63 gPARMS->SetDefaultParameter( "PI0FCALSKIM:MIN_R", MIN_R );
64 gPARMS->SetDefaultParameter( "PI0FCALSKIM:MAX_DT", MAX_DT );
65 gPARMS->SetDefaultParameter( "PI0FCALSKIM:MAX_ETOT", MAX_ETOT );
66 gPARMS->SetDefaultParameter( "PI0FCALSKIM:MIN_BLOCKS", MIN_BLOCKS );
67 gPARMS->SetDefaultParameter( "PI0FCALSKIM:WRITE_ROOT", WRITE_ROOT );
68 gPARMS->SetDefaultParameter( "PI0FCALSKIM:WRITE_EVIO", WRITE_EVIO );
69 */
70}
71
72//------------------
73// ~JEventProcessor_pi0fcalskim (Destructor)
74//------------------
75JEventProcessor_pi0fcalskim::~JEventProcessor_pi0fcalskim()
76{
77
78}
79
80//------------------
81// init
82//------------------
83jerror_t JEventProcessor_pi0fcalskim::init(void)
84{
85 dEventWriterEVIO = NULL__null;
86
87 num_epics_events = 0;
88/*
89 if( ! ( WRITE_ROOT || WRITE_EVIO ) ){
90
91 cerr << "No output mechanism has been specified." << endl;
92 return UNRECOVERABLE_ERROR;
93 }
94
95 if( WRITE_ROOT ){
96
97 japp->RootWriteLock();
98
99 m_tree = new TTree( "cluster", "Cluster Tree for Pi0 Calibration" );
100 m_tree->Branch( "nClus", &m_nClus, "nClus/I" );
101 m_tree->Branch( "hit0", m_hit0, "hit0[nClus]/I" );
102 m_tree->Branch( "px", m_px, "px[nClus]/F" );
103 m_tree->Branch( "py", m_py, "py[nClus]/F" );
104 m_tree->Branch( "pz", m_pz, "pz[nClus]/F" );
105
106 m_tree->Branch( "nHit", &m_nHit, "nHit/I" );
107 m_tree->Branch( "chan", m_chan, "chan[nHit]/I" );
108 m_tree->Branch( "e", m_e, "e[nHit]/F" );
109
110 japp->RootUnLock();
111 }
112*/
113 return NOERROR;
114}
115
116//------------------
117// brun
118//------------------
119jerror_t JEventProcessor_pi0fcalskim::brun(JEventLoop *eventLoop, int32_t runnumber)
120{
121 eventLoop->GetSingle(dEventWriterEVIO);
122
123 return NOERROR;
124}
125
126//------------------
127// evnt
128//------------------
129jerror_t JEventProcessor_pi0fcalskim::evnt(JEventLoop *loop, uint64_t eventnumber)
130{
131
132 vector< const DFCALShower* > locFCALShowers;
133 vector< const DVertex* > kinfitVertex;
134 loop->Get(locFCALShowers);
135 loop->Get(kinfitVertex);
136
137 vector< const DTrackTimeBased* > locTrackTimeBased;
138 loop->Get(locTrackTimeBased);
139
140 vector < const DFCALShower * > matchedShowers;
141
142 // always write out BOR events
143 if(loop->GetJEvent().GetStatusBit(kSTATUS_BOR_EVENT)) {
144 //jout << "Found BOR!" << endl;
145 dEventWriterEVIO->Write_EVIOEvent( loop, "pi0bcalskim" );
146 return NOERROR;
147 }
148
149 // write out the first few EPICS events to save run number & other meta info
150 if(loop->GetJEvent().GetStatusBit(kSTATUS_EPICS_EVENT) && (num_epics_events<5)) {
151 //jout << "Found EPICS!" << endl;
152 dEventWriterEVIO->Write_EVIOEvent( loop, "pi0bcalskim" );
153 num_epics_events++;
154 return NOERROR;
155 }
156
157
158 bool Candidate = false;
159
160 Double_t kinfitVertexX = 0.0, kinfitVertexY = 0.0, kinfitVertexZ = 0.0;
161
162 for (unsigned int i = 0 ; i < kinfitVertex.size(); i++)
163 {
164 kinfitVertexX = kinfitVertex[i]->dSpacetimeVertex.X();
165 kinfitVertexY = kinfitVertex[i]->dSpacetimeVertex.Y();
166 kinfitVertexZ = kinfitVertex[i]->dSpacetimeVertex.Z();
167 }
168
169 DVector3 norm(0.0,0.0,-1);
170 DVector3 pos,mom;
171 // Double_t radius = 0;
172 //japp->RootWriteLock();
173 Double_t p;
174 for (unsigned int i=0; i < locTrackTimeBased.size() ; ++i){
175 for (unsigned int j=0; j< locFCALShowers.size(); ++j){
176
177 Double_t x = locFCALShowers[j]->getPosition().X();
178 Double_t y = locFCALShowers[j]->getPosition().Y();
179 // Double_t z = locFCALShowers[j]->getPosition().Z() ;
180 //cout << "Z: " << z << endl;
181 //DVector3 pos_FCAL(x,y,625.406);
182 //for LH2 target
183 //DVector3 pos_FCAL(0,0,625.406);
184
185 DVector3 pos_FCAL(0,0,638);
186 //at the end of the start counter; use this fall for fall '15 data
187 // DVector3 pos_FCAL(0,0,692);
188 //DVector3 pos_FCAL(0.0,0.0,650);
189 //std::cout<<"i: "<< i<< " j: "<< j << " z: "<<z<< endl;
190 // if (locTrackTimeBased[i]->rt->GetIntersectionWithPlane(pos_FCAL,norm,pos,mom)==NOERROR)
191 if (locTrackTimeBased[i]->rt->GetIntersectionWithPlane(pos_FCAL,norm,pos,mom,NULL__null,NULL__null,NULL__null,SYS_FCAL)==NOERROR)
192 {
193 // Double_t dX= TMath::Abs(pos.X() - x);
194 // Double_t dY= TMath::Abs(pos.Y() - y);
195 // Double_t dZ= TMath::Abs(pos.Z() - z);
196 Double_t trkmass = locTrackTimeBased[i]->mass();
197 Double_t FOM = TMath::Prob(locTrackTimeBased[i]->chisq, locTrackTimeBased[i]->Ndof);
198 // radius = sqrt(pos.X()*pos.X() + pos.Y()*pos.Y());
199 // Double_t Eshwr = locFCALShowers[j]->getEnergy();
200 p = locTrackTimeBased[i]->momentum().Mag();
Value stored to 'p' is never read
201 // cout<<"p: "<<p<<endl;
202 // Double_t dZ = TMath::Abs(pos.Z() - z);
203 Double_t dRho = sqrt(((pos.X() - x)*(pos.X() - x)) + ((pos.Y() - y)* (pos.Y() - y)));
204 // std::cout<<"i: "<< i<< " j: "<< j << " dRho " <<dRho << endl;
205 //if(dX < 20 && dY < 20 && trkmass < 0.15 && dRho < 50 && FOM > 0.01) {
206 if(trkmass < 0.15 && dRho < 5 && FOM > 0.01 ) {
207 matchedShowers.push_back(locFCALShowers[j]);
208 // matchedTracks.push_back(locTrackTimeBased[i]);
209 // printf ("Matched event=%d, i=%d, j=%d, p=%f, Ztrk=%f Zshr=%f, Xtrk=%f, Xshr=%f, Ytrk=%f, Yshr=%f\n",locEventNumber,i,j,p,
210 // pos.Z(),z,pos.X(),x,pos.Y(),y);
211 // break;
212
213 }
214 }
215
216 }
217 }
218
219 for(unsigned int i=0; i<locFCALShowers.size(); i++)
220 {
221 if (find(matchedShowers.begin(), matchedShowers.end(),locFCALShowers[i]) != matchedShowers.end()) continue;
222
223 const DFCALShower *s1 = locFCALShowers[i];
224
225 vector<const DFCALCluster*> associated_clusters1;
226
227 s1->Get(associated_clusters1);
228 Double_t dx1 = s1->getPosition().X() - kinfitVertexX;
229 Double_t dy1 = s1->getPosition().Y() - kinfitVertexY;
230 Double_t dz1 = s1->getPosition().Z() - kinfitVertexZ;
231 Double_t R1 = sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
232 Double_t E1 = s1->getEnergy();
233 Double_t t1 = s1->getTime();
234 TLorentzVector sh1_p(E1*dx1/R1, E1*dy1/R1, E1*dz1/R1, E1);
235
236 for(unsigned int j=i+1; j<locFCALShowers.size(); j++)
237 {
238 const DFCALShower *s2 = locFCALShowers[j];
239 if (find(matchedShowers.begin(), matchedShowers.end(),s2) != matchedShowers.end()) continue;
240
241 vector<const DFCALCluster*> associated_clusters2;
242 s2->Get(associated_clusters2);
243 Double_t dx2 = s2->getPosition().X() - kinfitVertexX;
244 Double_t dy2 = s2->getPosition().Y() - kinfitVertexY;
245 Double_t dz2 = s2->getPosition().Z() - kinfitVertexZ;
246 Double_t R2 = sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2);
247 Double_t E2 = s2->getEnergy();
248 Double_t t2 = s2->getTime();
249
250 TLorentzVector sh2_p(E2*dx2/R2, E2*dy2/R2, E2*dz2/R2, E2);
251 TLorentzVector ptot = sh1_p+sh2_p;
252 Double_t inv_mass = ptot.M();
253
254 Candidate |= (E1 > 0.5 && E2 > 0.5 && s1->getPosition().Pt() > 20*k_cm && s2->getPosition().Pt() > 20*k_cm && (fabs (t1-t2) < 10) && (inv_mass<0.30) ) ;
255 }
256 }
257
258 if( Candidate ){
259
260 if( WRITE_EVIO ){
261
262 dEventWriterEVIO->Write_EVIOEvent( loop, "pi0fcalskim" );
263 }
264 }
265
266
267 /*
268 vector< const DFCALCluster* > clusterVec;
269 loop->Get( clusterVec );
270
271 if( clusterVec.size() < 2 ) return NOERROR;
272
273 bool hasCandidate = false;
274 double eTot = 0;
275
276 for( vector< const DFCALCluster*>::const_iterator clus1Itr = clusterVec.begin();
277 clus1Itr != clusterVec.end(); ++clus1Itr ){
278
279 eTot += (**clus1Itr).getEnergy();
280
281 for( vector< const DFCALCluster*>::const_iterator clus2Itr = clus1Itr + 1;
282 clus2Itr != clusterVec.end(); ++clus2Itr ){
283
284 const DFCALCluster& clusL =
285 ( (**clus1Itr).getEnergy() > (**clus2Itr).getEnergy() ?
286 (**clus2Itr) : (**clus1Itr) );
287
288 const DFCALCluster& clusH =
289 ( (**clus1Itr).getEnergy() > (**clus2Itr).getEnergy() ?
290 (**clus1Itr) : (**clus2Itr) );
291
292 double clusLX = clusL.getCentroid().X();
293 double clusLY = clusL.getCentroid().Y();
294 double rL = sqrt( clusLX * clusLX + clusLY * clusLY );
295 double eL = clusL.getEnergy();
296 double tL = clusL.getTime();
297 int nHitL = clusL.GetHits().size();
298
299 double clusHX = clusH.getCentroid().X();
300 double clusHY = clusH.getCentroid().Y();
301 double rH = sqrt( clusHX * clusHX + clusHY * clusHY );
302 double eH = clusH.getEnergy();
303 double tH = clusH.getTime();
304 int nHitH = clusH.GetHits().size();
305
306 DVector3 clusLMom = clusL.getCentroid();
307 clusLMom.SetMag( eL );
308
309 DVector3 clusHMom = clusH.getCentroid();
310 clusHMom.SetMag( eH );
311
312 double dt = fabs( tL - tH );
313
314 DLorentzVector gamL( clusLMom, clusLMom.Mag() );
315 DLorentzVector gamH( clusHMom, clusHMom.Mag() );
316
317 double mass = ( gamL + gamH ).M();
318
319 hasCandidate |=
320 ( ( eL > MIN_E ) &&
321 ( dt < MAX_DT ) &&
322 ( rL > MIN_R ) && ( rH > MIN_R ) &&
323 ( nHitL >= MIN_BLOCKS ) && ( nHitH >= MIN_BLOCKS ) &&
324 ( mass > MIN_MASS ) && ( mass < MAX_MASS ) );
325 }
326 }
327
328 if( hasCandidate && ( eTot < MAX_ETOT ) ){
329
330 if( WRITE_EVIO ){
331
332 dEventWriterEVIO->Write_EVIOEvent( loop, "pi0fcalskim" );
333 }
334
335 if( WRITE_ROOT ){
336
337 japp->RootWriteLock();
338 writeClustersToRoot( clusterVec );
339 japp->RootUnLock();
340 }
341 }
342*/
343 return NOERROR;
344}
345
346//------------------
347// erun
348//------------------
349jerror_t JEventProcessor_pi0fcalskim::erun(void)
350{
351 // This is called whenever the run number changes, before it is
352 // changed to give you a chance to clean up before processing
353 // events from the next run number.
354 return NOERROR;
355}
356
357//------------------
358// Fin
359//------------------
360jerror_t JEventProcessor_pi0fcalskim::fini(void)
361{
362 // Called before program exit after event processing is finished.
363 return NOERROR;
364}
365
366
367/*
368void
369JEventProcessor_pi0fcalskim::writeClustersToRoot( const vector< const DFCALCluster* > clusVec ){
370
371 // this code must run serially -- obtain a lock before
372 // entering this function
373
374 m_nHit = 0;
375 m_nClus = 0;
376
377 // hit and cluster indices
378 int& iH = m_nHit;
379 int& iC = m_nClus;
380
381 for( vector< const DFCALCluster* >::const_iterator clusItr = clusVec.begin();
382 clusItr != clusVec.end(); ++clusItr ){
383
384 // if we exceed max clusters abort writing for this event
385 if( iC == kMaxClus ) return;
386
387 const DFCALCluster& clus = (**clusItr);
388
389 if( ( clus.getCentroid().Perp() < 20*k_cm ) ||
390 ( clus.getEnergy() < 1*k_GeV ) ||
391 ( clus.GetHits().size() < 2 ) ) continue;
392
393 DVector3 gamMom = clus.getCentroid();
394 gamMom.SetMag( clus.getEnergy() );
395
396 m_hit0[iC] = iH;
397 m_px[iC] = gamMom.X();
398 m_py[iC] = gamMom.Y();
399 m_pz[iC] = gamMom.Z();
400
401 const vector<DFCALCluster::DFCALClusterHit_t>& hits = clus.GetHits();
402
403 for( unsigned int i = 0; i < hits.size(); ++i ){
404
405 // if we exceed max hits abort this event and return
406 if( iH == kMaxHits ) return;
407
408 m_chan[iH] = hits[i].ch;
409 m_e[iH] = hits[i].E;
410 ++iH;
411 }
412 ++iC;
413 }
414
415 m_tree->Fill();
416}
417*/