/[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.14 - (show annotations) (download)
Thu Sep 10 12:52:51 2009 UTC (15 years, 4 months ago) by mocchiut
Branch: MAIN
Changes since 1.13: +1 -0 lines
Debug printout added

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

  ViewVC Help
Powered by ViewVC 1.1.23