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

  ViewVC Help
Powered by ViewVC 1.1.23