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

  ViewVC Help
Powered by ViewVC 1.1.23