/[PAMELA software]/calo/flight/CaloPreSampler/src/CaloPreSampler.cpp
ViewVC logotype

Contents of /calo/flight/CaloPreSampler/src/CaloPreSampler.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.3 - (show annotations) (download)
Fri Aug 10 09:55:42 2007 UTC (17 years, 4 months ago) by mocchiut
Branch: MAIN
Changes since 1.2: +8 -0 lines
Small bugs fixed

1 /**
2 * \file CaloPreSampler.cpp
3 * \author Emiliano Mocchiutti (2007/07/18)
4 */
5 //
6 // headers
7 //
8 #include <CaloPreSampler.h>
9 //--------------------------------------
10 /**
11 * Default constructor
12 */
13 CaloPreSampler::CaloPreSampler(){
14 Clear();
15 };
16
17 CaloPreSampler::CaloPreSampler(PamLevel2 *l2p){
18 //
19 L2 = l2p;
20 //
21 if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
22 //
23 OBT = 0;
24 PKT = 0;
25 atime = 0;
26 //
27 // Default variables
28 //
29 event = new CaloLevel0();
30 cstrip = new CaloStrip(false);
31 // c1 = new CaloLevel1();
32 pcalo = new CaloLevel2();
33 N = 3;
34 debug = false;
35 sel = true;
36 cont = false;
37 emulate18 = true;
38 simulation = false;
39 //
40 Clear();
41 //
42 // loading magnetic field...
43 //
44 TrkLevel2 *trk = new TrkLevel2();
45 GL_PARAM *q4 = new GL_PARAM();
46 TSQLServer *dbc = 0;
47 TString host = "mysql://localhost/pamelaprod";
48 TString user = "anonymous";
49 TString psw = "";
50 const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
51 const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
52 const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
53 if ( !pamdbhost ) pamdbhost = "";
54 if ( !pamdbuser ) pamdbuser = "";
55 if ( !pamdbpsw ) pamdbpsw = "";
56 if ( strcmp(pamdbhost,"") ) host = pamdbhost;
57 if ( strcmp(pamdbuser,"") ) user = pamdbuser;
58 if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
59 dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
60 //
61 q4->Query_GL_PARAM(1,1,dbc);
62 printf(" Reading magnetic field maps at %s\n",(q4->PATH+q4->NAME).Data());
63 trk->LoadField(q4->PATH+q4->NAME);
64 //
65 };
66
67 void CaloPreSampler::Clear(){
68 //
69 pcalo->Clear();
70 //
71 };
72
73 void CaloPreSampler::Print(){
74 //
75 Process();
76 //
77 printf("========================================================================\n");
78 printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
79 printf(" debug [debug flag]:.. %i\n",debug);
80 printf(" simulation [simulation flag]:.. %i\n",simulation);
81 printf(" emulate18 [emulate dead plane 18]:.. %i\n",emulate18);
82 printf(" selection mode :.. %i\n",sel);
83 printf(" contamination mode :.. %i\n",cont);
84 printf(" pre-sampler planes :.. %i\n",N);
85 printf(" pcalo->qtot :.. %f\n",pcalo->qtot);
86 printf(" pcalo->nstrip :.. %i\n",pcalo->nstrip);
87 if ( pcalo->ntrk() > 0 ){
88 printf(" pcalo->track0->qtrack :.. %f\n",pcalo->GetCaloTrkVar(0)->qtrack);
89 printf(" pcalo->track0->dX0l :.. %f\n",pcalo->GetCaloTrkVar(0)->dX0l);
90 };
91 printf("========================================================================\n");
92 //
93 };
94
95 void CaloPreSampler::Delete(){
96 Clear();
97 delete pcalo;
98 //delete this;
99 };
100
101
102 void CaloPreSampler::Process(){
103 //
104 if ( !L2 ){
105 printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
106 printf(" ERROR: CaloPreSampler variables _NOT_ filled \n");
107 return;
108 };
109 //
110 Bool_t newentry = false;
111 //
112 if ( L2->IsORB() ){
113 if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || sel != ssel ){
114 newentry = true;
115 OBT = L2->GetOrbitalInfo()->OBT;
116 PKT = L2->GetOrbitalInfo()->pkt_num;
117 atime = L2->GetOrbitalInfo()->absTime;
118 ssel = sel;
119 };
120 } else {
121 newentry = true;
122 };
123 //
124 if ( !newentry ) return;
125 //
126 // Some variables
127 //
128 Int_t S3 = 0;
129 Int_t S2 = 0;
130 Int_t S12 = 0;
131 Int_t S11 = 0;
132 Float_t tmptrigty = -1.;
133 Bool_t trackanyway = true;
134 Float_t rigdefault = 50.;
135 Bool_t hZn = true;
136 Bool_t withtrk = true;
137 Bool_t st = true;
138 Int_t ntrkentry = 0;
139 TrkLevel2 *trk = L2->GetTrkLevel2();
140 Bool_t filled = false;
141 //
142 if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
143 //
144 this->Clear();
145 //
146 // find out if we have trkseqno = -1, -2 or -3
147 //
148 Bool_t m1 = false;
149 Bool_t m2 = false;
150 Bool_t m3 = false;
151 for (Int_t mm=0; mm < L2->GetCaloLevel2()->ntrk(); mm++ ){
152 if ( L2->GetCaloLevel2()->GetCaloTrkVar(mm)->trkseqno == -1 ) m1 = true;
153 if ( L2->GetCaloLevel2()->GetCaloTrkVar(mm)->trkseqno == -2 ) m2 = true;
154 if ( L2->GetCaloLevel2()->GetCaloTrkVar(mm)->trkseqno == -3 ) m3 = true;
155 };
156 //
157 if ( debug ) printf(" Fill estrip matrix needed to calculate variables \n");
158 //
159 // Fill the estrip matrix
160 //
161 memset(event->clevel1->estrip, 0, 2*22*96*sizeof(Float_t));
162 Int_t view = 0;
163 Int_t plane = 0;
164 Int_t strip = 0;
165 Float_t mip = 0.;
166 for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
167 //
168 mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
169 //
170 // Selection mode: fill the matrix only for plane < (22 - N) REMEMBER N = number of W planes to be used as presampler, ie if N = 2 then we want to use planes from 0 to 19
171 // included so plane < (22 - 2)
172 //
173 if ( sel ){
174 if ( plane < (22 - N) ){
175 //
176 if ( emulate18 && plane == (18 - N) ) mip = 0.;
177 event->clevel1->estrip[strip][plane][view] = mip;
178 //
179 };
180 };
181 //
182 // Contamination mode: fill the matrix only for planes from N to 22 but shift all planes up to the first one
183 //
184 if ( cont ){
185 if ( plane >= N ){
186 //
187 if ( emulate18 && plane == (18 + N) ) mip = 0.;
188 event->clevel1->estrip[strip][(plane-N)][view] = mip;
189 //
190 };
191 };
192 //
193 };
194 //
195 // if data comes from the simulation we must use mechanical alignment parameters (default is flight parameters)
196 //
197 if ( simulation ){
198 cstrip->UseMechanicalAlig();
199 };
200 //
201 // Set alignment parameter
202 //
203 event->clevel1->xalig = cstrip->GetXalig();
204 event->clevel1->yalig = cstrip->GetYalig();
205 event->clevel1->zalig = cstrip->GetZalig();
206 //
207 event->clevel1->emin = 0.7;
208 //
209 // in case of the contamination mode we must play with the Z alignment in order to have the correct track in the calo since we have moved the planes up...
210 //
211 if ( cont ){
212 if ( !(N%2) ){
213 event->clevel1->reverse = 0; // if the number of planes is even we have taken away a full module no need to do anything strange...
214 event->clevel1->zalig -= (N/2) * (8.09 + 10.09);
215 } else {
216 event->clevel1->reverse = 1; // if the number of planes is odd we have taken away half a module, we need to reverse silicon planes shifting
217 event->clevel1->zalig -= ((N+1)/2) * 8.09 + ((N-1)/2) * 10.09;
218 };
219 };
220 if ( debug ) printf(" xalig = %f \n",event->clevel1->xalig);
221 if ( debug ) printf(" yalig = %f \n",event->clevel1->yalig);
222 if ( debug ) printf(" zalig = %f \n",event->clevel1->zalig);
223 //
224 if ( debug ) printf(" Calculate variables as done in CaloCore, N = %i \n",N);
225 //
226 // Calculate variables
227 //
228 //
229 // use only N W planes
230 //
231 event->clevel1->npla = 22-N;
232 //
233 S3 = 0;
234 S2 = 0;
235 S12 = 0;
236 S11 = 0;
237 S3 = L2->GetTrigLevel2()->patterntrig[2];
238 S2 = L2->GetTrigLevel2()->patterntrig[3];
239 S12 = L2->GetTrigLevel2()->patterntrig[4];
240 S11 = L2->GetTrigLevel2()->patterntrig[5];
241 if ( L2->GetTrigLevel2()->patterntrig[1] & (1<<0) ) tmptrigty = 1.;
242 if ( L2->GetTrigLevel2()->patterntrig[0] ) tmptrigty = 2.;
243 if ( S3 || S2 || S12 || S11 ) tmptrigty = 0.;
244 if ( !(L2->GetTrigLevel2()->patterntrig[1] & (1<<0)) && !L2->GetTrigLevel2()->patterntrig[0] && !S3 && !S2 && !S12 && !S11 ) tmptrigty = 1.;
245 event->clevel2->trigty = tmptrigty;
246 //
247 // do we have at least one track from the tracker? this check has been disabled
248 //
249 event->clevel1->good2 = 1;
250 //
251 // copy variables calculated during calibration process which is skipped here...
252 //
253 event->clevel2->good = L2->GetCaloLevel2()->good;
254 memcpy(event->clevel2->perr,L2->GetCaloLevel2()->perr,sizeof(L2->GetCaloLevel2()->perr));
255 memcpy(event->clevel2->swerr,L2->GetCaloLevel2()->swerr,sizeof(L2->GetCaloLevel2()->swerr));
256 memcpy(event->clevel2->crc,L2->GetCaloLevel2()->crc,sizeof(L2->GetCaloLevel2()->crc));
257 event->clevel2->selftrigger = L2->GetCaloLevel2()->selftrigger;
258 //
259 // Calculate variables common to all tracks (qtot, nstrip, etc.)
260 //
261 if ( debug ) printf("1 Call GetCommonVar() \n");
262 event->GetCommonVar();
263 //
264 // Fill common variables
265 //
266 if ( debug ) printf("1 Call FillCommonVar() \n");
267 event->FillCommonVar(NULL,pcalo);
268 //
269 // Calculate variables related to tracks only if we have at least one track (from selftrigger and/or tracker)
270 //
271 ntrkentry = 0;
272 //
273 filled = false;
274 //
275 // Run over tracks (tracker or calorimeter )
276 //
277 if ( withtrk ){
278 //
279 for (Int_t nt=0; nt < trk->ntrk(); nt++){
280 //
281 event->clevel1->good2 = 1;
282 //
283 TrkTrack *ptt = trk->GetStoredTrack(nt);
284 //
285 event->clevel1->trkchi2 = 0;
286 //
287 // Copy the alpha vector in the input structure
288 //
289 for (Int_t e = 0; e < 5 ; e++){
290 event->clevel1->al_p[e][0] = ptt->al[e];
291 };
292 //
293 // Get tracker related variables for this track
294 //
295 if ( debug ) printf("track %i Call GetTrkVar() \n",nt);
296 event->GetTrkVar();
297 if ( debug ) printf(" event->clevel2->dX0l %f \n",event->clevel2->dX0l);
298 //
299 // Save tracker track sequence number
300 //
301 event->trkseqno = nt;
302 //
303 // Copy values in the class ca from the structure clevel2
304 //
305 if ( debug ) printf("track %i Call FillTrkVar() \n",nt);
306 event->FillTrkVar(pcalo,ntrkentry);
307
308
309 ntrkentry++;
310 filled = true;
311 //
312 }; // loop on all the tracks
313 };
314 //
315 // if no tracks found but there is the possibility to have a good track we should try to calculate anyway the track related variables using the calorimeter
316 // fit of the track (to be used for example when TRK is off due to any reason like IPM3/5 off).
317 // here we make an event selection so it must be done very carefully...
318 //
319 // conditions are: 0) no track from the tracker 1) we have a track fit both in x and y 2) no problems with calo for this event 3) no selftrigger event
320 //
321 // if ( trackanyway && !filled && event->clevel2->npcfit[0] >= 2 && event->clevel2->npcfit[1] >= 2 && event->clevel2->good != 0 && event->clevel2->trigty < 2. ){
322 if ( trackanyway && m3 ){
323 if ( debug ) printf(" Event with a track not fitted by the tracker \n");
324 //
325 // Disable "track mode" in the fortran routine
326 //
327 event->clevel1->good2 = 0;
328 event->clevel1->riginput = rigdefault;
329 if ( debug ) printf(" Using as default rigidity: %f \n",event->clevel1->riginput);
330 //
331 // We have a selftrigger event to analyze.
332 //
333 for (Int_t e = 0; e < 5 ; e++){
334 event->clevel1->al_p[e][0] = 0.;
335 event->clevel1->al_p[e][1] = 0.;
336 };
337 event->clevel1->trkchi2 = 0;
338 //
339 if ( debug ) printf("-3 a Call GetTrkVar() \n");
340 event->GetTrkVar();
341 //
342 // if we had no problem (clevel1->good2 = 0, NOTICE zero, not one in this mode!), fill and go on
343 //
344 if ( event->clevel1->good2 == 0 ) {
345 //
346 // In selftrigger mode the trkentry variable is set to -1
347 //
348 event->trkseqno = -3;
349 //
350 // Copy values in the class ca from the structure clevel2
351 //
352 if ( debug ) printf("-3 a Call FillTrkVar() \n");
353 event->FillTrkVar(pcalo,ntrkentry);
354 ntrkentry++;
355 filled = true;
356 //
357 } else {
358 if ( debug ) printf(" Selftrigger: problems with event \n");
359 };
360 //
361 };
362 //
363 // Call high energy nuclei routine
364 //
365 // if ( hZn && event->clevel2->trigty >= 2. ){
366 if ( hZn && m2 ){
367 if ( debug ) printf(" Calling selftrigger high energy nuclei routine \n");
368 //
369 // Disable "track mode" in the fortran routine
370 //
371 event->clevel1->good2 = 0;
372 //
373 // Set high energy nuclei flag to one
374 //
375 event->clevel1->hzn = 1;
376 event->clevel1->riginput = rigdefault;
377 //
378 // We have a selftrigger event to analyze.
379 //
380 for (Int_t e = 0; e < 5 ; e++){
381 event->clevel1->al_p[e][0] = 0.;
382 event->clevel1->al_p[e][1] = 0.;
383 };
384 event->clevel1->trkchi2 = 0;
385 //
386 if ( debug ) printf("-2 a Call GetTrkVar() \n");
387 event->GetTrkVar();
388 //
389 // if we had no problem (clevel1->good2 = 0, NOTICE zero, not one in this mode!), fill and go on
390 //
391 if ( event->clevel1->good2 == 0 ) {
392 //
393 // In selftrigger mode the trkentry variable is set to -1
394 //
395 event->trkseqno = -2;
396 //
397 // Copy values in the class ca from the structure clevel2
398 //
399 if ( debug ) printf("-2 a Call FillTrkVar() \n");
400 event->FillTrkVar(pcalo,ntrkentry);
401 ntrkentry++;
402 filled = true;
403 //
404 } else {
405 if ( debug ) printf(" Selftrigger: problems with event \n");
406 };
407 //
408 };
409 //
410 // self trigger event
411 //
412 // if ( st && event->clevel2->trigty >= 2. ){
413 if ( st && m1 ){
414 if ( debug ) printf(" Selftrigger event \n");
415 //
416 // Disable "track mode" in the fortran routine
417 //
418 event->clevel1->good2 = 0;
419 //
420 // disable high enery nuclei flag;
421 //
422 event->clevel1->hzn = 0;
423 //
424 // We have a selftrigger event to analyze.
425 //
426 for (Int_t e = 0; e < 5 ; e++){
427 event->clevel1->al_p[e][0] = 0.;
428 event->clevel1->al_p[e][1] = 0.;
429 };
430 event->clevel1->trkchi2 = 0;
431 //
432 if ( debug ) printf("-1 a Call GetTrkVar() \n");
433 event->GetTrkVar();
434 //
435 // if we had no problem (clevel2->good = 0, NOTICE zero, not one in selftrigger mode!), fill and go on
436 //
437 if ( event->clevel1->good2 == 0 ) {
438 //
439 // In selftrigger mode the trkentry variable is set to -1
440 //
441 event->trkseqno = -1;
442 //
443 // Copy values in the class ca from the structure clevel2
444 //
445 if ( debug ) printf("-1 a Call FillTrkVar() \n");
446 event->FillTrkVar(pcalo,ntrkentry);
447 ntrkentry++;
448 filled = true;
449 //
450 } else {
451 if ( debug ) printf(" Selftrigger: problems with event \n");
452 };
453 };
454 //
455 // Clear structures used to communicate with fortran
456 //
457 event->ClearStructs();
458 //
459 //
460 //
461 if ( debug ) this->Print();
462 if ( debug ) printf(" exit \n");
463 //
464 };

  ViewVC Help
Powered by ViewVC 1.1.23