/[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.12 - (show annotations) (download)
Tue Feb 24 09:18:05 2009 UTC (15 years, 9 months ago) by pam-fi
Branch: MAIN
Changes since 1.11: +56 -4 lines
New method AddCaloTrkVar(float* al, int trktag)

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

  ViewVC Help
Powered by ViewVC 1.1.23