Bug Summary

File:plugins/Utilities/pi0fcalskim/JEventProcessor_pi0fcalskim.cc
Location:line 185, 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 if( ! ( WRITE_ROOT || WRITE_EVIO ) ){
88
89 cerr << "No output mechanism has been specified." << endl;
90 return UNRECOVERABLE_ERROR;
91 }
92
93 if( WRITE_ROOT ){
94
95 japp->RootWriteLock();
96
97 m_tree = new TTree( "cluster", "Cluster Tree for Pi0 Calibration" );
98 m_tree->Branch( "nClus", &m_nClus, "nClus/I" );
99 m_tree->Branch( "hit0", m_hit0, "hit0[nClus]/I" );
100 m_tree->Branch( "px", m_px, "px[nClus]/F" );
101 m_tree->Branch( "py", m_py, "py[nClus]/F" );
102 m_tree->Branch( "pz", m_pz, "pz[nClus]/F" );
103
104 m_tree->Branch( "nHit", &m_nHit, "nHit/I" );
105 m_tree->Branch( "chan", m_chan, "chan[nHit]/I" );
106 m_tree->Branch( "e", m_e, "e[nHit]/F" );
107
108 japp->RootUnLock();
109 }
110*/
111 return NOERROR;
112}
113
114//------------------
115// brun
116//------------------
117jerror_t JEventProcessor_pi0fcalskim::brun(JEventLoop *eventLoop, int32_t runnumber)
118{
119 eventLoop->GetSingle(dEventWriterEVIO);
120
121 return NOERROR;
122}
123
124//------------------
125// evnt
126//------------------
127jerror_t JEventProcessor_pi0fcalskim::evnt(JEventLoop *loop, uint64_t eventnumber)
128{
129
130 vector< const DFCALShower* > locFCALShowers;
131 vector< const DVertex* > kinfitVertex;
132 loop->Get(locFCALShowers);
133 loop->Get(kinfitVertex);
134
135 vector< const DTrackTimeBased* > locTrackTimeBased;
136 loop->Get(locTrackTimeBased);
137
138 vector < const DFCALShower * > matchedShowers;
139
140
141
142
143 bool Candidate = false;
144
145 Double_t kinfitVertexX = 0.0, kinfitVertexY = 0.0, kinfitVertexZ = 0.0;
146
147 for (unsigned int i = 0 ; i < kinfitVertex.size(); i++)
148 {
149 kinfitVertexX = kinfitVertex[i]->dSpacetimeVertex.X();
150 kinfitVertexY = kinfitVertex[i]->dSpacetimeVertex.Y();
151 kinfitVertexZ = kinfitVertex[i]->dSpacetimeVertex.Z();
152 }
153
154 DVector3 norm(0.0,0.0,-1);
155 DVector3 pos,mom;
156 // Double_t radius = 0;
157 //japp->RootWriteLock();
158 Double_t p;
159 for (unsigned int i=0; i < locTrackTimeBased.size() ; ++i){
160 for (unsigned int j=0; j< locFCALShowers.size(); ++j){
161
162 Double_t x = locFCALShowers[j]->getPosition().X();
163 Double_t y = locFCALShowers[j]->getPosition().Y();
164 // Double_t z = locFCALShowers[j]->getPosition().Z() ;
165 //cout << "Z: " << z << endl;
166 //DVector3 pos_FCAL(x,y,625.406);
167 //for LH2 target
168 //DVector3 pos_FCAL(0,0,625.406);
169
170 DVector3 pos_FCAL(0,0,638);
171 //at the end of the start counter; use this fall for fall '15 data
172 // DVector3 pos_FCAL(0,0,692);
173 //DVector3 pos_FCAL(0.0,0.0,650);
174 //std::cout<<"i: "<< i<< " j: "<< j << " z: "<<z<< endl;
175 // if (locTrackTimeBased[i]->rt->GetIntersectionWithPlane(pos_FCAL,norm,pos,mom)==NOERROR)
176 if (locTrackTimeBased[i]->rt->GetIntersectionWithPlane(pos_FCAL,norm,pos,mom,NULL__null,NULL__null,NULL__null,SYS_FCAL)==NOERROR)
177 {
178 // Double_t dX= TMath::Abs(pos.X() - x);
179 // Double_t dY= TMath::Abs(pos.Y() - y);
180 // Double_t dZ= TMath::Abs(pos.Z() - z);
181 Double_t trkmass = locTrackTimeBased[i]->mass();
182 Double_t FOM = TMath::Prob(locTrackTimeBased[i]->chisq, locTrackTimeBased[i]->Ndof);
183 // radius = sqrt(pos.X()*pos.X() + pos.Y()*pos.Y());
184 // Double_t Eshwr = locFCALShowers[j]->getEnergy();
185 p = locTrackTimeBased[i]->momentum().Mag();
Value stored to 'p' is never read
186 // cout<<"p: "<<p<<endl;
187 // Double_t dZ = TMath::Abs(pos.Z() - z);
188 Double_t dRho = sqrt(((pos.X() - x)*(pos.X() - x)) + ((pos.Y() - y)* (pos.Y() - y)));
189 // std::cout<<"i: "<< i<< " j: "<< j << " dRho " <<dRho << endl;
190 //if(dX < 20 && dY < 20 && trkmass < 0.15 && dRho < 50 && FOM > 0.01) {
191 if(trkmass < 0.15 && dRho < 5 && FOM > 0.01 ) {
192 matchedShowers.push_back(locFCALShowers[j]);
193 // matchedTracks.push_back(locTrackTimeBased[i]);
194 // 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,
195 // pos.Z(),z,pos.X(),x,pos.Y(),y);
196 // break;
197
198 }
199 }
200
201 }
202 }
203
204 for(unsigned int i=0; i<locFCALShowers.size(); i++)
205 {
206 if (find(matchedShowers.begin(), matchedShowers.end(),locFCALShowers[i]) != matchedShowers.end()) continue;
207
208 const DFCALShower *s1 = locFCALShowers[i];
209
210 vector<const DFCALCluster*> associated_clusters1;
211
212 s1->Get(associated_clusters1);
213 Double_t dx1 = s1->getPosition().X() - kinfitVertexX;
214 Double_t dy1 = s1->getPosition().Y() - kinfitVertexY;
215 Double_t dz1 = s1->getPosition().Z() - kinfitVertexZ;
216 Double_t R1 = sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
217 Double_t E1 = s1->getEnergy();
218 Double_t t1 = s1->getTime();
219 TLorentzVector sh1_p(E1*dx1/R1, E1*dy1/R1, E1*dz1/R1, E1);
220
221 for(unsigned int j=i+1; j<locFCALShowers.size(); j++)
222 {
223 const DFCALShower *s2 = locFCALShowers[j];
224 if (find(matchedShowers.begin(), matchedShowers.end(),s2) != matchedShowers.end()) continue;
225
226 vector<const DFCALCluster*> associated_clusters2;
227 s2->Get(associated_clusters2);
228 Double_t dx2 = s2->getPosition().X() - kinfitVertexX;
229 Double_t dy2 = s2->getPosition().Y() - kinfitVertexY;
230 Double_t dz2 = s2->getPosition().Z() - kinfitVertexZ;
231 Double_t R2 = sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2);
232 Double_t E2 = s2->getEnergy();
233 Double_t t2 = s2->getTime();
234
235 TLorentzVector sh2_p(E2*dx2/R2, E2*dy2/R2, E2*dz2/R2, E2);
236 TLorentzVector ptot = sh1_p+sh2_p;
237 Double_t inv_mass = ptot.M();
238
239 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) ) ;
240 }
241 }
242
243 if( Candidate ){
244
245 if( WRITE_EVIO ){
246
247 dEventWriterEVIO->Write_EVIOEvent( loop, "pi0fcalskim" );
248 }
249 }
250
251
252 /*
253 vector< const DFCALCluster* > clusterVec;
254 loop->Get( clusterVec );
255
256 if( clusterVec.size() < 2 ) return NOERROR;
257
258 bool hasCandidate = false;
259 double eTot = 0;
260
261 for( vector< const DFCALCluster*>::const_iterator clus1Itr = clusterVec.begin();
262 clus1Itr != clusterVec.end(); ++clus1Itr ){
263
264 eTot += (**clus1Itr).getEnergy();
265
266 for( vector< const DFCALCluster*>::const_iterator clus2Itr = clus1Itr + 1;
267 clus2Itr != clusterVec.end(); ++clus2Itr ){
268
269 const DFCALCluster& clusL =
270 ( (**clus1Itr).getEnergy() > (**clus2Itr).getEnergy() ?
271 (**clus2Itr) : (**clus1Itr) );
272
273 const DFCALCluster& clusH =
274 ( (**clus1Itr).getEnergy() > (**clus2Itr).getEnergy() ?
275 (**clus1Itr) : (**clus2Itr) );
276
277 double clusLX = clusL.getCentroid().X();
278 double clusLY = clusL.getCentroid().Y();
279 double rL = sqrt( clusLX * clusLX + clusLY * clusLY );
280 double eL = clusL.getEnergy();
281 double tL = clusL.getTime();
282 int nHitL = clusL.GetHits().size();
283
284 double clusHX = clusH.getCentroid().X();
285 double clusHY = clusH.getCentroid().Y();
286 double rH = sqrt( clusHX * clusHX + clusHY * clusHY );
287 double eH = clusH.getEnergy();
288 double tH = clusH.getTime();
289 int nHitH = clusH.GetHits().size();
290
291 DVector3 clusLMom = clusL.getCentroid();
292 clusLMom.SetMag( eL );
293
294 DVector3 clusHMom = clusH.getCentroid();
295 clusHMom.SetMag( eH );
296
297 double dt = fabs( tL - tH );
298
299 DLorentzVector gamL( clusLMom, clusLMom.Mag() );
300 DLorentzVector gamH( clusHMom, clusHMom.Mag() );
301
302 double mass = ( gamL + gamH ).M();
303
304 hasCandidate |=
305 ( ( eL > MIN_E ) &&
306 ( dt < MAX_DT ) &&
307 ( rL > MIN_R ) && ( rH > MIN_R ) &&
308 ( nHitL >= MIN_BLOCKS ) && ( nHitH >= MIN_BLOCKS ) &&
309 ( mass > MIN_MASS ) && ( mass < MAX_MASS ) );
310 }
311 }
312
313 if( hasCandidate && ( eTot < MAX_ETOT ) ){
314
315 if( WRITE_EVIO ){
316
317 dEventWriterEVIO->Write_EVIOEvent( loop, "pi0fcalskim" );
318 }
319
320 if( WRITE_ROOT ){
321
322 japp->RootWriteLock();
323 writeClustersToRoot( clusterVec );
324 japp->RootUnLock();
325 }
326 }
327*/
328 return NOERROR;
329}
330
331//------------------
332// erun
333//------------------
334jerror_t JEventProcessor_pi0fcalskim::erun(void)
335{
336 // This is called whenever the run number changes, before it is
337 // changed to give you a chance to clean up before processing
338 // events from the next run number.
339 return NOERROR;
340}
341
342//------------------
343// Fin
344//------------------
345jerror_t JEventProcessor_pi0fcalskim::fini(void)
346{
347 // Called before program exit after event processing is finished.
348 return NOERROR;
349}
350
351
352/*
353void
354JEventProcessor_pi0fcalskim::writeClustersToRoot( const vector< const DFCALCluster* > clusVec ){
355
356 // this code must run serially -- obtain a lock before
357 // entering this function
358
359 m_nHit = 0;
360 m_nClus = 0;
361
362 // hit and cluster indices
363 int& iH = m_nHit;
364 int& iC = m_nClus;
365
366 for( vector< const DFCALCluster* >::const_iterator clusItr = clusVec.begin();
367 clusItr != clusVec.end(); ++clusItr ){
368
369 // if we exceed max clusters abort writing for this event
370 if( iC == kMaxClus ) return;
371
372 const DFCALCluster& clus = (**clusItr);
373
374 if( ( clus.getCentroid().Perp() < 20*k_cm ) ||
375 ( clus.getEnergy() < 1*k_GeV ) ||
376 ( clus.GetHits().size() < 2 ) ) continue;
377
378 DVector3 gamMom = clus.getCentroid();
379 gamMom.SetMag( clus.getEnergy() );
380
381 m_hit0[iC] = iH;
382 m_px[iC] = gamMom.X();
383 m_py[iC] = gamMom.Y();
384 m_pz[iC] = gamMom.Z();
385
386 const vector<DFCALCluster::DFCALClusterHit_t>& hits = clus.GetHits();
387
388 for( unsigned int i = 0; i < hits.size(); ++i ){
389
390 // if we exceed max hits abort this event and return
391 if( iH == kMaxHits ) return;
392
393 m_chan[iH] = hits[i].ch;
394 m_e[iH] = hits[i].E;
395 ++iH;
396 }
397 ++iC;
398 }
399
400 m_tree->Fill();
401}
402*/