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

  ViewVC Help
Powered by ViewVC 1.1.23