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

  ViewVC Help
Powered by ViewVC 1.1.23