This is a software test for reconstruction algorithm. The data used is 1000 events from monte carlo simulation. There are two sets of data. The physical values without the effect of detector transfer function and the data with the detector (the data will be acquired in the experiment). There is a problem in this especific test. In the experiment, when the particle decays, it is not considered in the reconstruction. In this case, there are combination of parents and offsprings.
1 // ------------------------------------------------------------------------
2 // recoevt.cc
3 // Description:
4 // Class recoevt - reconstruct parent particles from tracks and gammas
5 //
6 //------------------------------------------------------------------------
7 #include "BaBar/BaBar.hh"
8
9 //-----------------------
10 // This Class's Header --
11 //-----------------------
12 #include "BetaMiniUser/recoevt.hh"
13
14 //-------------
15 // C Headers --
16 //-------------
17 #include <assert.h>
18
19 //---------------
20 // C++ Headers --
21 //---------------
22 #include <iostream.h>
23 #include <math.h>
24
25 //-------------------------------
26 // Collaborating Class Headers --
27 //-------------------------------
28 #include "CLHEP/Alist/ConstAList.h" // produces a list
29 #include "CLHEP/Alist/ConstAIterator.h" // tool to iterate throeugh a list
30 #include "HepTuple/TupleManager.h"
31 #include "HepTuple/Histogram.h"
32 #include "HepTuple/Tuple.h"
33 #include "AbsEvent/AbsEvent.hh" //
34 #include "AbsEvent/getTmpAList.hh" // for making lists
35 #include "AbsEnv/AbsEnv.hh"
36 #include "GenEnv/GenEnv.hh"
37 #include "PDT/Pdt.hh" // needed for particle data table
38 #include "PDT/PdtPid.hh" // the same
39 #include "TrkBase/TrkRecoTrk.hh" // give you access to the track itself
40 // There you will find DC information
41 #include "TrkBase/TrkFit.hh"
42 #include "PDT/PdtEntry.hh" //
43 #include "BetaMicroAdapter/BtaCalQual.hh" // EMC information
44 #include "AbsParm/AbsParmIfdStrKey.hh"
45 #include "Beta/EventInfo.hh" // General event information
46 #include "Beta/BtaCandidate.hh" // Candidates are what everything is about
47 #include "BetaCoreTools/BtaMcAssoc.hh"
48 #include "ErrLogger/ErrLog.hh"
49
50 //-----------------------------------------------------------------------
51 // Local Macros, Typedefs, Structures, Unions and Forward Declarations --
52 //-----------------------------------------------------------------------
53
54 static const char rcsid[] = "$Id: recoevt.cc,v 1.7 2002/06/21 23:25:28 jtinslay Exp $";
55
56 //----------------
57 // Constructors --
58 //----------------
59 recoevt::recoevt(
60 const char* const theName,
61 const char* const theDescription )
62 : AppModule( theName, theDescription )
63 , _eventInfoList("eventInfoList",this,"Default") //we produce different
64 , _btaPionList("pionCandidates",this,"piDefault") //lists of candidates
65 , _btaCalorList("calorCandidates",this,"CalorNeutral") // We fill them with
66 , _btaTrackList("trackCandidates",this,"ChargedTracks")//default
67 , _btaGamList( "btaGamList", this, "CalorNeutral") // gamma
68 , _btaTruthList("truthCandidates",this,"MCTruth") //identification.
69 // This is handeled in myAnalysis.tcl
70 // You have to add the needed modules in AppUserBuild.cc too!
71 // (I have deleted muons from the candidate lists. You may try
72 // to make a muon list and loop over it)
73 {
74 commands()->append(& _eventInfoList);
75 commands()->append(& _btaPionList);
76 commands()->append(& _btaCalorList);
77 commands()->append(& _btaTrackList);
78 commands()->append(& _btaGamList);
79 commands()->append(& _btaTruthList);
80 }
81
82 //--------------
83 // Destructor --
84 //--------------
85 recoevt::~recoevt( ) {}
86
87 //--------------
88 // Operations --
89 //--------------
90
91 // The begin(AppJob*) member function is run before any events are
92 // processed. In this analysis, it opens the output histogram file
93 // and then books a number of histograms
94
95 AppResult
96 recoevt::beginJob( AbsEvent* anEvent )
97 {
98 _n_read_events = 0 ;
99
100 HepTupleManager* manager = gblEnv->getGen()->ntupleManager();
101 assert(manager != 0);
102
103 _massaMC2 = manager->histogram("Massa MC - comb 2", 12000, 0., 12. );
104 _massaSD2 = manager->histogram("Massa SD - comb 2", 12000, 0., 12. );
105 _massaMC3 = manager->histogram("Massa MC - comb 3", 12000, 0., 12. );
106 _massaSD3 = manager->histogram("Massa SD - comb 3", 12000, 0., 12. );
107 _massaMC4 = manager->histogram("Massa MC - comb 4", 12000, 0., 12. );
108 _massaSD4 = manager->histogram("Massa SD - comb 4", 12000, 0., 12. );
109 _massaMC5 = manager->histogram("Massa MC - comb 5", 12000, 0., 12. );
110 _massaSD5 = manager->histogram("Massa SD - comb 5", 12000, 0., 12. );
111 _massaMC6 = manager->histogram("Massa MC - comb 6", 12000, 0., 12. );
112 _massaSD6 = manager->histogram("Massa SD - comb 6", 12000, 0., 12. );
113
114 return AppResult::OK;
115 }
116
117 const int numpart=200;
118 int contpart,estado,lista[numpart];
119 float px[numpart],py[numpart],pz[numpart],energia[numpart];
120
121 void recoevt::comb(int init, int numtermos, int itera, int tam, int *lista)
122 {
123 int i,j,xlista[numpart];
124 float momx,momy,momz,ener,massa;
125
126 if((numtermos-itera)==1) {
127 for(i=lista[numtermos-2]+1;i<tam;i++) {
128 momx=0.0;
129 momy=0.0;
130 momz=0.0;
131 ener=0.0;
132 for(j=0;j<numtermos-1;j++) {
133 momx+=px[lista[j]];
134 momy+=py[lista[j]];
135 momz+=pz[lista[j]];
136 ener+=energia[lista[j]];
137 }
138 momx+=px[i];
139 momy+=py[i];
140 momz+=pz[i];
141 ener+=energia[i];
142 massa=sqrt(ener*ener-(momx*momx+momy*momy+momz*momz));
143 if(estado==0)
144 switch (numtermos) {
145 case 2:
146 _massaMC2->accumulate(massa);
147 break;
148 case 3:
149 _massaMC3->accumulate(massa);
150 break;
151 case 4:
152 _massaMC4->accumulate(massa);
153 break;
154 case 5:
155 _massaMC5->accumulate(massa);
156 break;
157 case 6:
158 _massaMC6->accumulate(massa);
159 break;
160 }
161 else
162 switch (numtermos) {
163 case 2:
164 _massaSD2->accumulate(massa);
165 break;
166 case 3:
167 _massaSD3->accumulate(massa);
168 break;
169 case 4:
170 _massaSD4->accumulate(massa);
171 break;
172 case 5:
173 _massaSD5->accumulate(massa);
174 break;
175 case 6:
176 _massaSD6->accumulate(massa);
177 break;
178 }
179 }
180 } else {
181 for(i=init;i<tam;i++) {
182 for(j=0;j<itera;j++)
183 xlista[j]=lista[j];
184 xlista[itera]=i;
185 comb(i+1,numtermos,itera+1,tam,xlista);
186 }
187 }
188 }
189
190 AppResult
191 recoevt::event( AbsEvent* anEvent ) {
192
193 int i;
194 _n_read_events++ ;
195 HepPoint mcPriVtx = HepPoint(0.,0.,0.);
196
197 // Get the event summary info object from the EventInfo list
198 HepAList<EventInfo>* eventList = Ifd<HepAList<EventInfo> >::get(anEvent,
199 _eventInfoList.value());
200 if (eventList == 0) {
201 ErrMsg(fatal) << "Event list is null. Something is wrong in the job configuration."
202 << endmsg;
203 }
204
205 // get list of charged tracks
206 HepAList<BtaCandidate>* trackList = Ifd<HepAList<BtaCandidate> >::get(anEvent,
207 _btaTrackList.value());
208 if (trackList == 0) {
209 ErrMsg(fatal) << "Charged track list is null. Something is wrong in the job configuration."
210 << endmsg;
211 }
212
213 // get list of MC truth particles
214 HepAList<BtaCandidate>* mcList = Ifd<HepAList<BtaCandidate> >::get(anEvent,
215 _btaTruthList.value());
216 if (mcList == 0) {
217 ErrMsg(fatal) << "MC list is null. Something is wrong in the job configuration."
218 << endmsg;
219 }
220
221 // get the default truth mapping
222 BtaMcAssoc* truthMap = Ifd< BtaMcAssoc >::get( anEvent, "GHit" );
223 assert(truthMap != 0); // should always be present, even if non-functional
224 // I don't know what's that.
225
226 // Loop over MC particles
227 HepConstAListIterator<BtaCandidate> iterMC(*mcList);
228 const BtaCandidate* mc_ptr;
229 int trknum=0;
230 while ( 0 != ( mc_ptr = iterMC() ) ) {
231
232 px[trknum]=mc_ptr->p4().x();
233 py[trknum]=mc_ptr->p4().y();
234 pz[trknum]=mc_ptr->p4().z();
235 energia[trknum]=mc_ptr->energy();
236
237 trknum++;
238 }
239
240 contpart=trknum;
241 lista[0]=0;
242 estado=0;
243 for(i=2;i<=4;i++)
244 comb(0,i,0,contpart,lista);
245
246 // Loop over track candidates
247 HepConstAListIterator<BtaCandidate> itertrack(*trackList);
248 const BtaCandidate* track_ptr;
249 contpart=0;
250 while ( 0 != ( track_ptr = itertrack() ) ) {
251
252 px[contpart]=track_ptr->p4().x();
253 py[contpart]=track_ptr->p4().y();
254 pz[contpart]=track_ptr->p4().z();
255 energia[contpart]=track_ptr->energy();
256
257 contpart++;
258 } // end: while ( 0 != ( track_ptr = itertrack() ) )
259
260
261 // photon list
262
263 HepAList<BtaCandidate>* gamList =
264 Ifd<HepAList<BtaCandidate> >::get(anEvent, _btaGamList.value() );
265
266 HepAListIterator<BtaCandidate> gamIter(*gamList);
267 BtaCandidate* theGamma;
268 while (theGamma = gamIter() ) {
269
270 px[contpart]=theGamma->p4().x();
271 py[contpart]=theGamma->p4().y();
272 pz[contpart]=theGamma->p4().z();
273 energia[contpart]=theGamma->energy();
274
275 contpart++;
276 } // end: while ( 0 != ( theGamma = gamIter() ) )
277
278 lista[0]=0;
279 estado=1;
280 for(i=2;i<=4;i++)
281 comb(0,i,0,contpart,lista);
282
283 return AppResult::OK;
284 }
285
286 AppResult
287 recoevt::endJob( AbsEvent* anEvent )
288 {
289 return AppResult::OK;
290 }
1 //--------------------------------------------------------------------------
2 // recoevt.hh
3 //
4 // Description:
5 // Class recoevt - reconstruct the parent particle from traks and gammas
6 //
7 //------------------------------------------------------------------------
8
9 #ifndef RECOEVT_HH
10 #define RECOEVT_HH
11
12 //----------------------
13 // Base Class Headers --
14 //----------------------
15 #include "Framework/APPModule.hh"
16 #include "Framework/AbsParmBool.hh"
17 #include "AbsParm/AbsParmIfdStrKey.hh"
18
19 //------------------------------------
20 // Collaborating Class Declarations --
21 //------------------------------------
22 class HepTupleManager;
23 class HepHistogram;
24 class HepTuple;
25
26 // ---------------------
27 // -- Class Interface --
28 // ---------------------
29
30 class recoevt : public AppModule {
31
32 //--------------------
33 // Instance Members --
34 //--------------------
35
36 public:
37
38 // Constructors
39 recoevt( const char* const theName, const char* const theDescription );
40
41 // Destructor
42 virtual ~recoevt( );
43
44 // Operations
45
46 virtual AppResult beginJob( AbsEvent* anEvent);
47 virtual AppResult event( AbsEvent* anEvent );
48 virtual AppResult endJob ( AbsEvent* anEvent );
49 virtual void comb(int init, int numtermos, int itera, int tam, int *lista);
50
51 protected:
52
53 AbsParmIfdStrKey _eventInfoList;
54 AbsParmIfdStrKey _btaCalorList;
55 AbsParmIfdStrKey _btaTrackList;
56 AbsParmIfdStrKey _btaGamList;
57 AbsParmIfdStrKey _btaTruthList;
58 AbsParmIfdStrKey _btaPionList;
59
60 private:
61 int _n_read_events; // this is our event counter
62
63 HepHistogram* _massaMC2;
64 HepHistogram* _massaSD2;
65 HepHistogram* _massaMC3;
66 HepHistogram* _massaSD3;
67 HepHistogram* _massaMC4;
68 HepHistogram* _massaSD4;
69 HepHistogram* _massaMC5;
70 HepHistogram* _massaSD5;
71 HepHistogram* _massaMC6;
72 HepHistogram* _massaSD6;
73
74 };
75
76 #endif
1 #------------------------------------------------------------------------------
2 # MomEnH.tcl
3 #------------------------------------------------------------------------------
4 # always source the error logger early in your main tcl script
5 sourceFoundFile ErrLogger/ErrLog.tcl
6 sourceFoundFile FrameScripts/FwkCfgVar.tcl
7 sourceFoundFile FrameScripts/talkto.tcl
8 # Disable the use of envvars
9 set ProdTclOnly true
10
11 # set the error logging level to 'warning'. If you encounter a configuration
12 # error you can get more information using 'trace'
13 ErrLoggingLevel warning
14
15 #
16 # Turn off some specialty items
17 #
18 module disable NTrk
19 module disable MomEnH
20 module disable PEntp
21 module disable colPEntp
22 module disable EnPBip
23 module disable LisEvt
24 module disable verMCTruth
25
26 ## allowed values of BetaMiniReadPersistence are (currently) "Kan", "Bdb"
27 ##
28 set BetaMiniReadPersistence Kan
29
30 ## allowed (non-expert) values of levelOfDetail are "micro", "cache", "extend"
31 ## or "refit"
32 ##
33 FwkCfgVar levelOfDetail "cache"
34
35 ## allowed values of ConfigPatch are "Run1", "Run2" or "MC". This MUST be set
36 ## consistent ## with your input data type or you will get INCONSISTENT OR
37 ## INCORRECT RESULTS
38 ##
39 FwkCfgVar ConfigPatch "MC"
40
41 ##
42 ## You can enter input collections two ways: either append them to a list, or
43 ## explicitly enter them in the input module. Do one or the other, BUT NOT
44 ## BOTH.
45 ## If inputList is set before executing btaMini.tcl, that will automatically
46 ## add the collections to the appropriate input module, otherwise make sure you
47 ## talk to the right one.
48 ##
49 ## lappend inputList collection1 collection2 ...
50 ##
51 ## OR THE FOLLOWING (choose the correct one based on persistence)
52 ##
53 ## talkto BdbEventInput {
54 ## talkto KanEventInput {
55 ## input add collection1
56 ## input add collection2
57 ## ...
58 ## }
59
60
61 ##
62 ## Set the number of events to run. If this isn't set, all events in the
63 ## input collections will be processed.
64 ##
65 FwkCfgVar NEvent
66
67 ## choose the flavor of ntuple to write (hbook or root) and the file name
68 ##
69 FwkCfgVar BetaMiniTuple "root"
70 FwkCfgVar histFileName "enphist.roo"
71
72 ## create Everything path and add core sequences to it. btaMiniPhysics is the
73 ## same as btaMini, just appending a few standard list generating modules. For
74 ## reading data with stored composites, you may have a conflict running
75 ## btaMiniPhyscs.tcl
76 ##
77 ## You can also run (most of) the PhysProdSequence, complete with its 3 gamma
78 ## conversion finders, etc. Consider disabling the portion of this sequence
79 ## that you do not need to save yourself some time. The BetaLumiSequence
80 ## and TagProd sequences are left off, as they otherwise cause problems.
81 ##
82 sourceFoundFile BetaMiniUser/btaMini.tcl
83 #sourceFoundFile BetaMiniUser/btaMiniPhysics.tcl
84 #sourceFoundFile BetaMiniUser/btaMiniPhysProdSequence.tcl
85 ## Add Analysis module
86 ##
87 path append Everything recoevt
88
89 ##
90 ## If your job has a tag-level filter, here is how you should run it
91 ## so as to avoid wasting time reading the mini when the tag filter fails
92 ## Here's a simple example that restricts to just multi-hadron events
93 ## on Kan input
94
95 module clone TagFilterByName TagBGFMultiHadron
96 module talk TagBGFMultiHadron
97 andList set BGFMultiHadron
98 assertIfMissing set true
99 exit
100
101 module talk EvtCounter
102 printFreq set 100000
103 exit
104
105 #sequence append BetaMiniReadSequence -a KanEventUpdateTag TagBGFMultiHadron
106
107 path list
108 if [info exists NEvent] {
109 ev begin -nev $NEvent
110 } else {
111 ev begin
112 }
113
114 ErrMsg trace "completed OK"
115 exit
Most of the peaks and resonances disappear from detector data. They are not reconstruct from the final stables particles because all neutrinos do not interact with the detector, and they carry a considerable amount of energy.
|
|
|
Feedback to: jamwer@hep.man.ac.uk |