File: | plugins/Utilities/pi0fcalskim/JEventProcessor_pi0fcalskim.cc |
Location: | line 200, column 4 |
Description: | Value stored to 'p' is never read |
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" |
11 | using 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 | |
33 | extern "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 | //------------------ |
44 | JEventProcessor_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 | //------------------ |
75 | JEventProcessor_pi0fcalskim::~JEventProcessor_pi0fcalskim() |
76 | { |
77 | |
78 | } |
79 | |
80 | //------------------ |
81 | // init |
82 | //------------------ |
83 | jerror_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 | //------------------ |
119 | jerror_t JEventProcessor_pi0fcalskim::brun(JEventLoop *eventLoop, int32_t runnumber) |
120 | { |
121 | eventLoop->GetSingle(dEventWriterEVIO); |
122 | |
123 | return NOERROR; |
124 | } |
125 | |
126 | //------------------ |
127 | // evnt |
128 | //------------------ |
129 | jerror_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 | //------------------ |
349 | jerror_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 | //------------------ |
360 | jerror_t JEventProcessor_pi0fcalskim::fini(void) |
361 | { |
362 | // Called before program exit after event processing is finished. |
363 | return NOERROR; |
364 | } |
365 | |
366 | |
367 | /* |
368 | void |
369 | JEventProcessor_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 | */ |