/[PAMELA software]/calo/flight/CaloPreSampler/src/CaloPreSampler.cpp
ViewVC logotype

Annotation of /calo/flight/CaloPreSampler/src/CaloPreSampler.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.7 - (hide annotations) (download)
Mon Dec 3 13:57:11 2007 UTC (17 years ago) by mocchiut
Branch: MAIN
Changes since 1.6: +2 -1 lines
Small bug in the latest feature fixed

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

  ViewVC Help
Powered by ViewVC 1.1.23