/[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.16 - (show annotations) (download)
Mon Dec 14 14:31:55 2009 UTC (15 years ago) by mocchiut
Branch: MAIN
Changes since 1.15: +3 -0 lines
Do not use plane 18x by default

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

  ViewVC Help
Powered by ViewVC 1.1.23