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

  ViewVC Help
Powered by ViewVC 1.1.23