/[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.10 - (hide annotations) (download)
Wed May 28 08:31:46 2008 UTC (16 years, 6 months ago) by mocchiut
Branch: MAIN
Changes since 1.9: +2 -0 lines
Mask plane 18X itself when bool emulate18 is true

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

  ViewVC Help
Powered by ViewVC 1.1.23