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

  ViewVC Help
Powered by ViewVC 1.1.23