/[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.4 - (hide annotations) (download)
Fri Aug 31 15:36:35 2007 UTC (17 years, 3 months ago) by mocchiut
Branch: MAIN
Changes since 1.3: +5 -2 lines
Added standalone capabilities

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

  ViewVC Help
Powered by ViewVC 1.1.23