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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1.18 - (hide annotations) (download)
Thu Dec 18 10:31:00 2008 UTC (16 years, 10 months ago) by mocchiut
Branch: MAIN
Changes since 1.17: +4 -4 lines 
Re-apply ROOT 5.20 compilation bug fixed
1 mocchiut 1.1 #include <CaloNuclei.h>
2 malvezzi 1.14 #include <TGraph.h>
3     #include <TSpline.h>
4     #include <TMVA/TSpline1.h>
5 mocchiut 1.1
6     //--------------------------------------
7     /**
8     * Default constructor
9     */
10     CaloNuclei::CaloNuclei(){
11     Clear();
12     };
13    
14     CaloNuclei::CaloNuclei(PamLevel2 *l2p){
15     //
16     Clear();
17     //
18     L2 = l2p;
19     //
20     if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
21     //
22     OBT = 0;
23     PKT = 0;
24     atime = 0;
25     N = 5;
26 mocchiut 1.2 R = 3;
27 mocchiut 1.1 //
28 mocchiut 1.4 debug = false;
29 malvezzi 1.14 // debug = true;
30 mocchiut 1.5 usetrack = true;
31 mocchiut 1.4 //
32 mocchiut 1.1 };
33    
34     void CaloNuclei::Clear(){
35     //
36 mocchiut 1.11 UN = 0;
37 mocchiut 1.3 tr = 0;
38 mocchiut 1.5 sntr = 0;
39 mocchiut 1.1 interplane = 0;
40     preq = 0.;
41     postq = 0.;
42 mocchiut 1.10 stdedx1 = 0.;
43     ethr = 0.;
44 mocchiut 1.1 dedx1 = 0.;
45     dedx3 = 0.;
46     qpremean = 0.;
47     qpremeanN = 0.;
48 malvezzi 1.14 maxrel = 0;
49     qNmin1 = 0;
50     qNmin1_w = 0;
51     charge_siegen1 = 0;
52     ZCalo_maxrel_b = 0;
53     ZCalo_dedx_b = 0;
54     ZCalo_dedx_defl= 0;
55     ZCalo_Nmin1_defl= 0;
56 mocchiut 1.1 //
57     multhit = false;
58     gap = false;
59     //
60     };
61    
62     void CaloNuclei::Print(){
63     //
64 mocchiut 1.2 Process();
65     //
66 mocchiut 1.3 printf("========================================================================\n");
67 mocchiut 1.5 printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack);
68 malvezzi 1.14 printf(" interplane [number of available dE/dx before interaction]:....... %i\n",interplane);
69     printf(" ethr [threshold used to determine interplane]:................... %f \n",ethr);
70     printf(" dedx1 [dE/dx from the first calorimeter plane]:.................. %f \n",dedx1);
71     printf(" stdedx1 [dE/dx from the first calorimeter plane standalone]:..... %f \n",stdedx1);
72     printf(" dedx3 [dE/dx (average) if the first 3 Si planes]:................ %f \n",dedx3);
73     printf(" multhit [true if interplane determined by multiple hits]:........ %i \n",multhit);
74     printf(" gap [true if interplane determined by a gap]:.................... %i \n",gap);
75     printf(" preq [total energy in MIP before the interaction plane]:......... %f \n",preq);
76     printf(" postq [total energy in MIP after the interaction plane]:......... %f \n",postq);
77     printf(" qpremean [truncated mean using 3 planes and 3 strips]:........... %f \n",qpremean);
78     printf(" N [no of used plane]:............................................ %i \n",N);
79     printf(" R [no strip used per plane ]:.................................... %i \n",R);
80     printf(" qpremeanN [truncated mean using N planes and R strips]:.......... %f \n",qpremeanN);
81     printf(" qNmin1 [truncated mean using N-1 planes and R strips]: .......... %f \n",qNmin1);
82     printf(" maxrel [dE/dx of strip with maximum release (I plane)]:.......... %f \n",maxrel);
83     printf(" ZCalo_maxrel_b [Z from maximum release in I Calo plane vs beta].. %f \n",ZCalo_maxrel_b);
84     printf(" ZCalo_dedx_b [Z from dedx in I Calo plane vs beta].. ............ %f \n",ZCalo_dedx_b);
85     printf(" ZCalo_dedx_defl [Z from dedx in I Calo plane vs deflection....... %f \n",ZCalo_dedx_defl);
86     printf(" ZCalo_Nmin1_defl [Z from truncated mean (N-1 pl) vs deflection].. %f \n",ZCalo_Nmin1_defl);
87 mocchiut 1.3 printf("========================================================================\n");
88 mocchiut 1.1 //
89     };
90    
91     void CaloNuclei::Delete(){
92     Clear();
93     //delete this;
94     };
95    
96    
97     void CaloNuclei::Process(){
98 mocchiut 1.3 Process(0);
99     };
100    
101     void CaloNuclei::Process(Int_t ntr){
102     //
103 mocchiut 1.1 if ( !L2 ){
104     printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
105     printf(" ERROR: CaloNuclei variables not filled \n");
106     return;
107     };
108     //
109     Bool_t newentry = false;
110     //
111     if ( L2->IsORB() ){
112 mocchiut 1.5 if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || ntr != sntr ){
113 mocchiut 1.1 newentry = true;
114     OBT = L2->GetOrbitalInfo()->OBT;
115     PKT = L2->GetOrbitalInfo()->pkt_num;
116     atime = L2->GetOrbitalInfo()->absTime;
117 mocchiut 1.5 sntr = ntr;
118 mocchiut 1.1 };
119     } else {
120     newentry = true;
121     };
122     //
123     if ( !newentry ) return;
124     //
125 mocchiut 1.3 tr = ntr;
126     //
127 mocchiut 1.4 if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
128 mocchiut 1.1 //
129     Clear();
130 mocchiut 1.11 //
131 mocchiut 1.6 if ( debug ) printf(" Always calculate stdedx1 \n");
132 mocchiut 1.1 //
133 malvezzi 1.14 // Always calculate stdedx1 and maxrel
134 mocchiut 1.1 //
135 malvezzi 1.14 Int_t cont=0;
136 mocchiut 1.1 Int_t view = 0;
137     Int_t plane = 0;
138     Int_t strip = 0;
139 mocchiut 1.5 Int_t indx = 0;
140     Float_t vfpl[96];
141     Int_t stfpl[96];
142     memset(vfpl, 0, 96*sizeof(Float_t));
143     memset(stfpl, 0, 96*sizeof(Int_t));
144 mocchiut 1.1 Float_t mip = 0.;
145 mocchiut 1.5 for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
146     //
147     mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
148     //
149     // put in vfpl vector the energy release on the first plane
150     //
151     if ( strip != -1 && view == 1 && plane == 0 ) {
152     stfpl[indx] = strip;
153     vfpl[indx] = mip;
154     indx++;
155     };
156     //
157     };
158     //
159 mocchiut 1.6 if ( debug ) printf(" find energy released along the strip of maximum on the first plane and on the two neighbour strips \n");
160     //
161 mocchiut 1.5 // find energy released along the strip of maximum on the first plane and on the two neighbour strips
162     //
163     if ( indx > 0 ){
164 mocchiut 1.12 Int_t mindx = (Int_t)TMath::LocMax(indx,vfpl);
165 mocchiut 1.5 for (Int_t ii=0; ii<indx; ii++){
166     if ( stfpl[ii] == stfpl[mindx] ) stdedx1 += vfpl[ii];
167 mocchiut 1.13 if ( (mindx-1)>=0 && stfpl[ii] == (stfpl[mindx]-1) ) stdedx1 += vfpl[ii];
168     if ( (mindx+1)<96 && stfpl[ii] == (stfpl[mindx]+1) ) stdedx1 += vfpl[ii];
169     // if ( (mindx-1)>=0 && stfpl[ii] == stfpl[mindx-1] ) stdedx1 += vfpl[ii];
170     // if ( (mindx+1)<96 && stfpl[ii] == stfpl[mindx+1] ) stdedx1 += vfpl[ii];
171 mocchiut 1.5 };
172 malvezzi 1.14 maxrel = vfpl[mindx];
173 mocchiut 1.5 } else {
174     stdedx1 = 0.;
175 malvezzi 1.14 maxrel = 0.;
176 mocchiut 1.5 };
177 malvezzi 1.14 // cout<<stdedx1<<" "<<maxrel<<"\n";
178 mocchiut 1.5 //
179 mocchiut 1.6 if ( debug ) printf(" if ( !usetrack ) return: usetrack %i ntr %i \n",usetrack,ntr);
180 mocchiut 1.5 //
181     //
182 mocchiut 1.8 // if ( !usetrack ) return;
183 mocchiut 1.5 //
184     PamTrack *ptrack = 0;
185     CaloTrkVar *track = 0;
186 mocchiut 1.8 //
187     if ( usetrack ){
188     if ( ntr >= 0 ){
189     ptrack = L2->GetTrack(ntr);
190     if ( ptrack ) track = ptrack->GetCaloTrack();
191     } else {
192     track = L2->GetCaloStoredTrack(ntr);
193     };
194     //
195     if ( !track && ntr >= 0 ){
196     printf(" ERROR: cannot find any track!\n");
197     printf(" ERROR: CaloNuclei variables not completely filled \n");
198     return;
199     };
200 mocchiut 1.5 } else {
201 mocchiut 1.8 if ( ntr >= 0 ){
202 mocchiut 1.9 if ( debug ) printf(" ERROR: you asked not to use a track but you are looking for track number %i !\n",ntr);
203     if ( debug ) printf(" ERROR: CaloNuclei variables not completely filled \n");
204 mocchiut 1.8 return;
205     };
206 mocchiut 1.5 };
207     //
208 mocchiut 1.1 // Float_t defethr = 6. * 0.90;
209     Float_t defethr = 6.25; // = (sqrt(9) - 0.5) ** 2.;
210     //
211     // Calculate dedx1 and dedx3
212     //
213     for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
214     //
215     mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
216     //
217 mocchiut 1.5 if ( ntr >= 0 ){
218     //
219     if ( strip != -1 &&
220     view == 1 &&
221     plane == 0 &&
222     ( strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) || strip == (track->tibar[0][1]) )
223     && true ){
224     dedx1 += mip;
225     };
226     if ( strip != -1 &&
227     (( view == 1 && ( plane == 0 || plane == 1 ) ) ||
228     ( view == 0 && plane == 0 )) &&
229     (( view == 0 && ( strip == track->tibar[0][0] || strip == (track->tibar[0][0]-1) || strip == (track->tibar[0][0]-2) )) ||
230     ( view == 1 && ( strip == track->tibar[0][1] || strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) )) ||
231     ( view == 1 && ( strip == track->tibar[1][1] || strip == (track->tibar[1][1]-1) || strip == (track->tibar[1][1]-2) ))) &&
232     true ){
233     dedx3 += mip;
234     };
235     } else {
236     //
237     if ( strip != -1 &&
238     view == 1 &&
239     plane == 0 &&
240     ( strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) || strip == (L2->GetCaloLevel2()->cibar[0][1]) )
241     && true ){
242     dedx1 += mip;
243     };
244     if ( strip != -1 &&
245     (( view == 1 && ( plane == 0 || plane == 1 ) ) ||
246     ( view == 0 && plane == 0 )) &&
247     (( view == 0 && ( strip == L2->GetCaloLevel2()->cibar[0][0] || strip == (L2->GetCaloLevel2()->cibar[0][0]-1) || strip == (L2->GetCaloLevel2()->cibar[0][0]-2) )) ||
248     ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[0][1] || strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) )) ||
249     ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[1][1] || strip == (L2->GetCaloLevel2()->cibar[1][1]-1) || strip == (L2->GetCaloLevel2()->cibar[1][1]-2) ))) &&
250     true ){
251     dedx3 += mip;
252     };
253 mocchiut 1.1 };
254     //
255     };
256     //
257     dedx3 /= 3.;
258     // Float_t mesethr = dedx1 * 0.90;
259     Float_t mesethr = 0.;
260     if ( dedx1 > 0. ) mesethr = (sqrt(dedx1) - 0.50)*(sqrt(dedx1) - 0.50);
261     Bool_t aldone = false;
262     //
263     retry:
264     //
265 mocchiut 1.4 if ( debug ) printf("retry\n");
266 mocchiut 1.1 //
267     interplane = 0;
268     //
269     ethr = TMath::Max(defethr,mesethr);
270     //
271     // Find the interaction plane "interplane"
272     //
273     Int_t gapth = 3;
274     Int_t nhit[2] = {0,0};
275     Int_t splane[2] = {-1,-1};
276     Int_t sview[2] = {-1,-1};
277     Int_t interpl[2] = {-1,-1};
278     Int_t interv[2] = {-1,-1};
279     Bool_t wmulthit[2] = {false,false};
280     Bool_t wgap[2] = {false,false};
281     Int_t ii = 0;
282     while ( ii<L2->GetCaloLevel1()->istrip ){
283     //
284     mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);
285     //
286 mocchiut 1.5 if ( ntr >= 0 ){
287     if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
288     ( strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) )
289     && true ){
290     if ( debug ) printf(" inside loop: ii %i mip %f view %i plane %i strip %i tibar %i nhit %i splane %i sview %i \n",ii,mip,view,plane,strip,track->tibar[plane][view]-1,nhit[view],splane[view],sview[view]);
291     interpl[view] = plane;
292     interv[view] = view;
293     if ( splane[view] != plane || sview[view] != view ){
294     if ( nhit[view] > 1 ){
295     wmulthit[view] = true;
296     // if ( splane[view] == -1 ) splane[view] = 0; //
297     // if ( sview[view] == -1 ) sview[view] = view; //
298     interpl[view] = splane[view];
299     interv[view] = sview[view];
300     };
301     if ( plane > splane[view]+gapth ){
302     wgap[view] = true;
303     // if ( splane[view] == -1 ) splane[view] = 0;//
304     // if ( sview[view] == -1 ) sview[view] = view; //
305     interpl[view] = splane[view];
306     interv[view] = sview[view];
307     };
308     splane[view] = plane;
309     sview[view] = view;
310     nhit[view] = 1;
311     } else {
312     nhit[view]++;
313     };
314     };
315     } else {
316     if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
317     ( strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) )
318     && true ){
319     if ( debug ) printf(" inside loop: ii %i mip %f view %i plane %i strip %i cibar %i nhit %i splane %i sview %i \n",ii,mip,view,plane,strip,L2->GetCaloLevel2()->cibar[plane][view]-1,nhit[view],splane[view],sview[view]);
320     interpl[view] = plane;
321     interv[view] = view;
322     if ( splane[view] != plane || sview[view] != view ){
323     if ( nhit[view] > 1 ){
324     wmulthit[view] = true;
325     // if ( splane[view] == -1 ) splane[view] = 0; //
326     // if ( sview[view] == -1 ) sview[view] = view; //
327     interpl[view] = splane[view];
328     interv[view] = sview[view];
329 mocchiut 1.1 };
330 mocchiut 1.5 if ( plane > splane[view]+gapth ){
331     wgap[view] = true;
332     // if ( splane[view] == -1 ) splane[view] = 0;//
333     // if ( sview[view] == -1 ) sview[view] = view; //
334     interpl[view] = splane[view];
335     interv[view] = sview[view];
336     };
337     splane[view] = plane;
338     sview[view] = view;
339     nhit[view] = 1;
340     } else {
341     nhit[view]++;
342 mocchiut 1.1 };
343     };
344     };
345     //
346     ii++;
347     //
348     };
349     //
350 mocchiut 1.4 if (debug ) printf("conversion interpl %i interv %i multhit %i interplane %i \n",interpl[0],interv[0],multhit,interplane);
351 mocchiut 1.1 Int_t winterplane[2] = {-1,-1};
352     //
353     for ( Int_t view = 0; view < 2; view++){
354     //
355     if ( nhit[view] > 1 && !wmulthit[view] && !wgap[view] ){
356     wmulthit[view] = true;
357     interpl[view] = splane[view];
358     interv[view] = sview[view];
359     };
360     //
361     if ( wmulthit[view] ) multhit = true;
362     if ( wgap[view] ) gap = true;
363     //
364     // convert view and plane number of interaction plane into number of available dE/dx measurements before the interaction plane
365     //
366     if ( interpl[view] >= 0 ) {
367     if ( interv[view] == 0 ){
368     winterplane[view] = (1 + interpl[view]) * 2;
369     } else {
370     winterplane[view] = (1 + interpl[view]) + (1 + interpl[view] - 1);
371     };
372     if ( wmulthit[view] ) winterplane[view]--;
373     };
374     };
375     if ( winterplane[0] > 0 && winterplane[1] > 0 ){
376     if ( multhit ){
377     interplane = TMath::Min(winterplane[0],winterplane[1]);
378     } else {
379     interplane = TMath::Max(winterplane[0],winterplane[1]);
380     };
381     } else {
382     if ( !winterplane[0] || !winterplane[1] ){
383     interplane = 0;
384     } else {
385     interplane = TMath::Max(winterplane[0],winterplane[1]);
386     };
387     };
388     //
389 mocchiut 1.4 if ( debug ) printf("2conversion interpl %i interv %i multhit %i interplane %i \n",interpl[1],interv[1],multhit,interplane);
390     if ( debug ) printf("3conversion winterpl0 %i winterpl1 %i \n",winterplane[0],winterplane[1]);
391 mocchiut 1.1 //
392     Int_t ipl = 0;
393     if ( interplane > 0 ){
394     //
395     // Calculate preq, postq, qpremean
396     //
397 malvezzi 1.14 cont++;
398 mocchiut 1.1 ii = 0;
399     Int_t ind = -1;
400     Int_t qsplane = -1;
401     Int_t qsview = -1;
402 mocchiut 1.2 Int_t ind2 = -1;
403     Int_t qsplane2 = -1;
404     Int_t qsview2 = -1;
405 mocchiut 1.1 Float_t qme[200];
406     memset(qme,0,200*sizeof(Float_t));
407 mocchiut 1.2 Float_t qme2[2112];
408     memset(qme2,0,2112*sizeof(Float_t));
409 mocchiut 1.1 //
410     while ( ii<L2->GetCaloLevel1()->istrip ){
411     //
412     mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);
413     //
414     if ( strip != -1 ){
415     if ( view == 0 ){
416     ipl = (1 + plane) * 2;
417     } else {
418     ipl = (1 + plane) + (1 + plane - 1 );
419     };
420     if ( ipl > interplane ){
421     postq += mip;
422     } else {
423     preq += mip;
424 mocchiut 1.5 if ( ntr >= 0 ){
425     if ( strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) ){
426     if ( qsplane != plane || qsview != view ){
427     qsplane = plane;
428     qsview = view;
429     ind++;
430     if ( debug && ind > 199 ) printf(" AAAGH!! \n");
431     qme[ind] = 0.;
432     };
433     qme[ind] += mip;
434 mocchiut 1.1 };
435 mocchiut 1.5 for ( Int_t ns = 0; ns < R ; ns++){
436     Int_t ms = track->tibar[plane][view] - 1 - ns + (R - 1)/2;
437     if ( strip == ms ){
438     if ( qsplane2 != plane || qsview2 != view ){
439     qsplane2 = plane;
440     qsview2 = view;
441     ind2++;
442     if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
443     qme2[ind2] = 0.;
444     };
445     qme2[ind2] += mip;
446     };
447     };
448     } else {
449     if ( strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) ){
450     if ( qsplane != plane || qsview != view ){
451     qsplane = plane;
452     qsview = view;
453     ind++;
454     if ( debug && ind > 199 ) printf(" AAAGH!! \n");
455     qme[ind] = 0.;
456 mocchiut 1.2 };
457 mocchiut 1.5 qme[ind] += mip;
458 mocchiut 1.2 };
459 mocchiut 1.5 for ( Int_t ns = 0; ns < R ; ns++){
460     Int_t ms = L2->GetCaloLevel2()->cibar[plane][view] - 1 - ns + (R - 1)/2;
461     if ( strip == ms ){
462     if ( qsplane2 != plane || qsview2 != view ){
463     qsplane2 = plane;
464     qsview2 = view;
465     ind2++;
466     if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
467     qme2[ind2] = 0.;
468     };
469     qme2[ind2] += mip;
470     };
471     };
472     };
473 mocchiut 1.1 };
474     //
475     };
476     //
477     ii++;
478     //
479     };
480 malvezzi 1.14
481    
482     //
483 mocchiut 1.1 // here we must calculate qpremean, order vector qme, select 3 lowest measurements and caculate the mean...
484     //
485 mocchiut 1.4 if ( debug ){
486     for (Int_t l=0; l < interplane; l++){
487     printf(" qme[%i] = %f \n",l,qme[l]);
488     };
489     };
490     //
491 mocchiut 1.1 Long64_t work[200];
492     ind = 0;
493     Int_t l = 0;
494 mocchiut 1.5 Int_t RN = 0;
495 mocchiut 1.1 Float_t qm = 0.;
496 mocchiut 1.2 Float_t qm2 = 0.;
497     //
498     Float_t qmt = ethr*0.8; // *0.9
499     //
500 mocchiut 1.1 Int_t uplim = TMath::Max(3,N);
501 malvezzi 1.14 Int_t uplim2 = interplane-1;
502 mocchiut 1.2 //
503 mocchiut 1.1 while ( l < uplim && ind < interplane ){
504 mocchiut 1.18 qm = TMath::KOrdStat((Long64_t)interplane,qme,(Long64_t)ind,work);
505 mocchiut 1.2 if ( qm >= qmt ){
506 mocchiut 1.5 if ( l < 3 ){
507     qpremean += qm;
508     RN++;
509     };
510 mocchiut 1.1 l++;
511 mocchiut 1.4 if ( debug ) printf(" value no %i qm %f qmt %f \n",l,qm,qmt);
512 mocchiut 1.2 };
513     ind++;
514     };
515     //
516 mocchiut 1.5 qpremean /= (Float_t)RN;
517 mocchiut 1.2 ind = 0;
518     l = 0;
519 mocchiut 1.5 RN = 0;
520 mocchiut 1.2 while ( l < uplim && ind < interplane ){
521 mocchiut 1.18 qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
522 mocchiut 1.2 if ( qm2 >= qmt ){
523 mocchiut 1.5 if ( l < N ){
524     qpremeanN += qm2;
525     RN++;
526     };
527 mocchiut 1.2 l++;
528 mocchiut 1.5 if ( debug ) printf(" qm2 value no %i qm %f qmt %f RN %i \n",l,qm2,qmt,RN);
529 mocchiut 1.1 };
530     ind++;
531     };
532 malvezzi 1.14 ////////////////////////////////////
533     //to calculate qNmin1///////////////
534     ///////////////////////////////////
535     //values under threshold
536     qm2=0;
537     ind = 0;
538     l = 0;
539     RN = 0;
540     S2=0;
541     while ( l < uplim2 && ind<interplane){
542 mocchiut 1.18 qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
543 malvezzi 1.14 if ( qm2 < qmt ) S2++;
544     ind++;
545     }
546     qm2=0;
547     ind = 0;
548     l = 0;
549     RN = 0;
550     while ( l < uplim2 && ind < interplane ){
551 mocchiut 1.18 qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
552 malvezzi 1.14 if ( qm2 >= qmt ){
553     if ( l < (interplane - 1 - S2)){
554     qNmin1_w += qm2;
555     RN++;
556     };
557     l++;
558     if ( debug ) printf(" qm2 value no %i qm %f qmt %f RN %i \n",l,qm2,qmt,RN);
559     };
560     ind++;
561     };
562 mocchiut 1.5 qpremeanN /= (Float_t)RN;
563 malvezzi 1.14 qNmin1_w /= (Float_t)RN;
564 mocchiut 1.11 UN = RN;
565 malvezzi 1.14 ///////set qNmin1 definition///////////
566     if (interplane==1 || interplane==2){
567     if (dedx1>0) qNmin1=dedx1;
568     else if (stdedx1>0) qNmin1=stdedx1;
569     }
570     else if (interplane > 2){
571     qNmin1 = qNmin1_w;
572     }
573     ////////////////////////////////////
574     //////////////////////////////////
575 mocchiut 1.1 //
576 mocchiut 1.4 if ( debug ) printf(" charge is %f \n",sqrt(qpremean));
577 mocchiut 1.1 //
578     if ( mesethr != ethr && interplane >= 3 && !aldone ){
579     Float_t mesethr2 = (sqrt(qpremean) - 0.50)*(sqrt(qpremean) - 0.50);
580     if ( mesethr2 < mesethr*0.90 ){
581     mesethr = (sqrt(dedx1) - 0.25)*(sqrt(dedx1) - 0.25);
582     } else {
583     mesethr = mesethr2;
584     };
585     aldone = true;
586 mocchiut 1.7 if ( mesethr > defethr ){
587     interplane = 0;
588     preq = 0.;
589     postq = 0.;
590     qpremean = 0.;
591     qpremeanN = 0.;
592 malvezzi 1.14 qNmin1 = 0;
593 mocchiut 1.7 multhit = false;
594     gap = false;
595     goto retry;
596     };
597 mocchiut 1.1 };
598     };
599 malvezzi 1.14
600    
601    
602     //=======================================================================
603     //=========== charge determination stdedx1 vs. beta ===============
604     //====================== Siegen method ===========================
605     //=======================================================================
606    
607     // Data from file Calo_Bands_New_7.dat
608     Float_t C0[9] = {0 , 1 , 2 , 3 , 4 , 5 , 6 , 8 , 90 };
609     Float_t B0[9] = {0 , -2.03769 , 7.61781 , 19.7098 , 60.5598 , 57.9226 , 14.8368 , -1358.83 , 8200 };
610     Float_t B1[9] = {0 , 0.0211274 , 9.32057e-010 , 4.47241e-07 , 1.44826e-06 , 2.6189e-05 , 0.00278178 , 55.5445 , 0 };
611     Float_t B2[9] = {0 , -3.91742 , -20.0359 , -16.3043 , -16.9471 , -14.4622 , -10.9594 , -2.38014 , 0 };
612     Float_t B3[9] = {0 , 11.1469 , -6.63105 , -27.8834 , -132.044 , -55.341 , 173.25 , 4115 , 0 };
613     Float_t B4[9] = {0 , -14.3465 , -0.485215 , 18.8122 , 117.533 , -14.0898 , -325.269 , -4388.89 , 0 };
614     Float_t B5[9] = {0 , 6.24281 , 3.96018 , 0 , -26.1881 , 42.9731 , 182.697 , 1661.01 , 0 };
615    
616     Float_t x1[9],y1[9];
617     Int_t n1 = 9;
618    
619     Float_t charge = 1000.;
620     Float_t beta = 100.;
621    
622     //------- First try track dependent beta
623     if( L2->GetTrkLevel2()->GetNTracks()>=1 ){
624     PamTrack *TRKtrack = L2->GetTrack(0);
625     if (fabs(TRKtrack->GetToFTrack()->beta[12]) < 100.) beta = fabs(TRKtrack->GetToFTrack()->beta[12]);
626     }
627     //------- If no beta found, try standalone beta
628     if (beta == 100.) {
629     ToFTrkVar *ttrack = L2->GetToFStoredTrack(-1);
630     beta = fabs(ttrack->beta[12]);
631     }
632    
633     if (beta<2.) { // it makes no sense to allow beta=5 or so...
634    
635     Float_t mip=0;
636     mip=stdedx1 ;
637    
638     if (mip>0) {
639    
640     Float_t betahelp = pow(beta, 1.8);
641     Float_t ym = mip*betahelp;
642     Float_t xb = beta;
643    
644     for ( Int_t jj=0; jj<9; jj++ ){
645     x1[jj] = B0[jj]+B1[jj]*pow(xb,B2[jj])+B3[jj]*xb+B4[jj]*xb*xb+B5[jj]*xb*xb*xb;
646     y1[jj] = C0[jj]*C0[jj] ;
647     }
648    
649     TGraph *gr1 = new TGraph(n1,x1,y1);
650     TSpline3 *spl1 = new TSpline3("grs",gr1); // use a cubic spline
651     Float_t chelp = spl1->Eval(ym);
652     charge = TMath::Sqrt(chelp);
653     gr1->Delete();
654     spl1->Delete();
655    
656     } // if (mip1>0)
657     } // beta < 100
658    
659    
660     charge_siegen1 = charge;
661    
662     //======================= end charge Siegen ===========================
663    
664    
665     // //=======================================================================
666     // //=========== charge determination Maximum release vs. beta ===============
667     // //====================== Rome method ===========================
668 marcelli 1.17 // //=======================================================================
669    
670     Float_t D0[9] = {0, 1, 2, 3 , 4 , 5 , 6, 8, 90};
671     Float_t E1[9] = {0, 500, 500, 923.553 , 659.842, 1113.97, 3037.25, 3034.84, 0};
672     Float_t E2[9] = {0, 11.0, 7.5, 6.92574 , 5.08865, 5.29349, 6.41442, 5.52969, 0};
673     Float_t E3[9] = {0, 1.2, 4, 9.7227 , 13.18, 23.5444, 38.2057, 63.6784, 80000};
674 malvezzi 1.14
675 marcelli 1.17 Float_t xx1[9],yy1[9];
676     n1 = 9;
677 malvezzi 1.14
678     charge = 1000.;
679     mip=0;
680    
681    
682     if (beta<2.) { // it makes no sense to allow beta=5 or so...
683    
684    
685     mip=maxrel;
686    
687     if (mip>0) {
688     Float_t ym = mip;
689     Float_t xb = beta;
690    
691     for ( Int_t jj=0; jj<n1; jj++ ){
692     xx1[jj] = E1[jj]*exp(-E2[jj]*xb)+E3[jj];
693     yy1[jj] = D0[jj]*D0[jj] ;
694     }
695    
696     TGraph *gr1 = new TGraph(n1,xx1,yy1);
697     TSpline3 *spl1 = new TSpline3("grs",gr1); // use a cubic spline
698     Float_t chelp = spl1->Eval(ym);
699     charge = TMath::Sqrt(chelp);
700     gr1->Delete();
701     spl1->Delete();
702    
703    
704     } // if (mip1>0)
705     } // beta < 100
706    
707    
708     ZCalo_maxrel_b = charge;
709    
710     //======================= end charge Rome: maxril vs beta ===========================
711    
712    
713    
714     // =======================================================================
715     // =========== charge determination dedx vs. beta ===============
716     // ====================== Rome method ===========================
717 marcelli 1.17 // =======================================================================
718    
719     Float_t F0[9] = {0., 1., 2., 3. ,4., 5. , 6., 8, 90};
720     Float_t G1[9] = {0, 500, 500, 642.935 , 848.684, 1346.05, 3238.82, 3468.6, 0};
721     Float_t G2[9] = {0, 11, 7.5, 6.2038 , 5.51723, 5.65265, 6.54089, 5.72723, 0};
722     Float_t G3[9] = {0, 1.2, 4, 9.2421 , 13.9858, 25.3912, 39.6332, 64.5674, 80000};
723 malvezzi 1.14
724    
725     charge = 1000.;
726     mip=0;
727    
728    
729     if (beta<2.) { // it makes no sense to allow beta=5 or so...
730    
731    
732     if( L2->GetTrkLevel2()->GetNTracks()>=1 ){
733     mip=dedx1;
734     }
735     if (mip==0) mip=stdedx1;
736    
737    
738     if (mip>0) {
739    
740     Float_t ym = mip;
741     Float_t xb = beta;
742    
743     for ( Int_t jj=0; jj<n1; jj++ ){
744     xx1[jj] = G1[jj]*exp(-G2[jj]*xb)+G3[jj];
745     yy1[jj] = F0[jj]*F0[jj] ;
746     }
747    
748     TGraph *gr1 = new TGraph(n1,xx1,yy1);
749     TSpline3 *spl1 = new TSpline3("grs",gr1); // use a cubic spline
750     Float_t chelp = spl1->Eval(ym);
751     charge = TMath::Sqrt(chelp);
752     gr1->Delete();
753     spl1->Delete();
754    
755     } //if (mip1>0)
756     } //beta < 100
757    
758     ZCalo_dedx_b = charge;
759    
760     //======================= end charge Rome: dedx vs beta ===========================
761    
762    
763     //=======================================================================
764     //=========== charge determination dedx vs. defl ===============
765     //====================== Rome method ===========================
766 marcelli 1.17 //=======================================================================
767 malvezzi 1.14
768 marcelli 1.17 Float_t H0[9] = {0, 1, 2, 3 , 4 , 5 , 6, 8, 90 };
769     Float_t I1[9] = {0, 3.5, 40, 56.1019, 101.673, 109.225, 150.599, 388.531, 0};
770     Float_t I2[9] = {0, -1, -13.6, -12.5581, -22.5543, -15.9823, -28.2207, -93.6871, 0};
771     Float_t I3[9] = {0, 1, 5.3, 11.6218, 19.664, 32.1817, 45.7527, 84.5992, 80000};
772 malvezzi 1.14
773    
774     charge = 1000.;
775     mip=0;
776     Float_t defl=0;
777    
778    
779     if (beta<2.) { // it makes no sense to allow beta=5 or so...
780    
781     if( L2->GetTrkLevel2()->GetNTracks()>=1 ){
782     PamTrack *TRKtrack = L2->GetTrack(0);
783     mip=dedx1;
784     if (mip==0) mip=stdedx1;
785     defl=TRKtrack->GetTrkTrack()->al[4];
786    
787    
788     if (mip>0 && defl<0.7 && defl>0) {
789    
790     Float_t ym = mip;
791     Float_t xb = defl;
792    
793     for ( Int_t jj=0; jj<n1; jj++ ){
794     xx1[jj] = I1[jj]*xb*xb+I2[jj]*xb+I3[jj];
795     yy1[jj] = H0[jj]*H0[jj] ;
796     }
797    
798     TGraph *gr1 = new TGraph(n1,xx1,yy1);
799     TSpline3 *spl1 = new TSpline3("grs",gr1); // use a cubic spline
800     Float_t chelp = spl1->Eval(ym);
801     charge = TMath::Sqrt(chelp);
802     gr1->Delete();
803     spl1->Delete();
804    
805     } // if (mip1>0 && defl<0.5 && defl>0)
806     }//Ntrack>=1
807     } //beta < 100
808    
809     ZCalo_dedx_defl = charge;
810    
811     //======================= end charge Rome: dedx vs defl ===========================
812    
813    
814     //============================================================================================
815     //=========== charge determination Truncated mean (N-1 planes) vs. defl ===================
816     //================================ Rome method ========================================
817 marcelli 1.17 //============================================================================================
818    
819     Float_t L0[9] = {0, 1, 2, 3 , 4 , 5 , 6, 8, 90};
820     Float_t M1[9] = {0, 3.5, 27, 63.0145, 120.504, 173.663, 245.33, 236.517, 0};
821     Float_t M2[9] = {0, -1, -10.6, -15.005, -31.0635, -39.4988, -60.5011, -46.3992, 0};
822     Float_t M3[9] = {0, 1, 7, 12.5037, 22.8652, 35.2907, 51.4678, 86.4155, 80000};
823 malvezzi 1.14
824     charge = 1000.;
825     mip=0;
826    
827    
828     if (beta<2.) { // it makes no sense to allow beta=5 or so...
829    
830     if( L2->GetTrkLevel2()->GetNTracks()>=1 ){
831     mip=qNmin1;
832    
833     if (mip>0 && defl<0.7 && defl>0) {
834    
835     Float_t ym = mip;
836     Float_t xb = defl;
837    
838     for ( Int_t jj=0; jj<n1; jj++ ){
839     xx1[jj] = M1[jj]*xb*xb+M2[jj]*xb+M3[jj];
840     yy1[jj] = L0[jj]*L0[jj] ;
841     }
842    
843     TGraph *gr1 = new TGraph(n1,xx1,yy1);
844     TSpline3 *spl1 = new TSpline3("grs",gr1); // use a cubic spline
845     Float_t chelp = spl1->Eval(ym);
846     charge = TMath::Sqrt(chelp);
847     gr1->Delete();
848     spl1->Delete();
849    
850     } // if (mip1>0 && defl<0.5 && defl>0)
851     }//Ntrack>=1
852     } //beta < 100
853    
854     ZCalo_Nmin1_defl = charge;
855    
856     //======================= end charge Rome: Nmin1 vs defl ===========================
857    
858    
859    
860    
861    
862 mocchiut 1.1 //
863 mocchiut 1.5 if ( debug ) this->Print();
864 mocchiut 1.4 if ( debug ) printf(" esci \n");
865 mocchiut 1.1 //
866     };
  ViewVC Help
Powered by ViewVC 1.1.23