File: | plugins/Utilities/pi0fcalskim/JEventProcessor_pi0fcalskim.cc |
Location: | line 185, 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 | 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 | //------------------ |
117 | jerror_t JEventProcessor_pi0fcalskim::brun(JEventLoop *eventLoop, int32_t runnumber) |
118 | { |
119 | eventLoop->GetSingle(dEventWriterEVIO); |
120 | |
121 | return NOERROR; |
122 | } |
123 | |
124 | //------------------ |
125 | // evnt |
126 | //------------------ |
127 | jerror_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 | //------------------ |
334 | jerror_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 | //------------------ |
345 | jerror_t JEventProcessor_pi0fcalskim::fini(void) |
346 | { |
347 | // Called before program exit after event processing is finished. |
348 | return NOERROR; |
349 | } |
350 | |
351 | |
352 | /* |
353 | void |
354 | JEventProcessor_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 | */ |