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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1.1.1 by mocchiut, Tue Apr 3 15:55:16 2007 UTC revision 1.11 by mocchiut, Mon Nov 26 08:48:24 2007 UTC
# Line 20  CaloNuclei::CaloNuclei(PamLevel2 *l2p){ Line 20  CaloNuclei::CaloNuclei(PamLevel2 *l2p){
20    PKT = 0;    PKT = 0;
21    atime = 0;    atime = 0;
22    N = 5;    N = 5;
23      R = 3;
24      //
25      debug = false;
26      usetrack = true;
27    //    //
28  };  };
29    
30  void CaloNuclei::Clear(){  void CaloNuclei::Clear(){
31    //    //
32      UN = 0;
33      tr = 0;
34      sntr = 0;
35    interplane = 0;    interplane = 0;
36    preq = 0.;    preq = 0.;
37    postq = 0.;    postq = 0.;
38      stdedx1 = 0.;
39      ethr = 0.;
40    dedx1 = 0.;    dedx1 = 0.;
41    dedx3 = 0.;    dedx3 = 0.;
42    qpremean = 0.;    qpremean = 0.;
# Line 40  void CaloNuclei::Clear(){ Line 49  void CaloNuclei::Clear(){
49    
50  void CaloNuclei::Print(){  void CaloNuclei::Print(){
51    //    //
52    printf("==============================\n");    Process();
53    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);    //
54    printf(" interplane: %i\n",interplane);    printf("========================================================================\n");
55    printf(" ethr: %f \n",ethr);    printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack);
56    printf(" dedx1: %f \n",dedx1);    printf(" interplane [number of available dE/dx before interaction]:.. %i\n",interplane);
57    printf(" dedx3: %f \n",dedx3);    printf(" ethr [threshold used to determine interplane]:.............. %f \n",ethr);
58    printf(" multhit: %i \n",multhit);    printf(" dedx1 [dE/dx from the first calorimeter plane]:............. %f \n",dedx1);
59    printf(" gap: %i \n",gap);    printf(" stdedx1 [dE/dx from the first calorimeter plane standalone]: %f \n",stdedx1);
60    printf(" preq: %f \n",preq);    printf(" dedx3 [dE/dx (average) if the first 3 Si planes]:........... %f \n",dedx3);
61    printf(" postq: %f \n",postq);    printf(" multhit [true if interplane determined by multiple hits]:... %i \n",multhit);
62    printf(" qpremean: %f \n",qpremean);    printf(" gap [true if interplane determined by a gap]:............... %i \n",gap);
63    printf(" N: %i \n",N);    printf(" preq [total energy in MIP before the interaction plane]:.... %f \n",preq);
64    printf(" qpremeanN: %f \n",qpremeanN);    printf(" postq [total energy in MIP after the interaction plane]:.... %f \n",postq);
65    printf("==============================\n");    printf(" qpremean [truncated mean using 3 planes and 3 strips]:...... %f \n",qpremean);
66      printf(" N [no of used plane]:....................................... %i \n",N);
67      printf(" R [no strip used per plane ]:............................... %i \n",R);
68      printf(" qpremeanN [truncated mean using N planes and R strips]:..... %f \n",qpremeanN);
69      printf("========================================================================\n");
70    //    //
71  };  };
72    
# Line 64  void CaloNuclei::Delete(){ Line 77  void CaloNuclei::Delete(){
77    
78    
79  void CaloNuclei::Process(){  void CaloNuclei::Process(){
80    //    Process(0);
81    };
82    
83    void CaloNuclei::Process(Int_t ntr){
84      //  
85    if ( !L2 ){    if ( !L2 ){
86      printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");      printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
87      printf(" ERROR: CaloNuclei variables not filled \n");      printf(" ERROR: CaloNuclei variables not filled \n");
# Line 74  void CaloNuclei::Process(){ Line 91  void CaloNuclei::Process(){
91    Bool_t newentry = false;    Bool_t newentry = false;
92    //    //
93    if ( L2->IsORB() ){    if ( L2->IsORB() ){
94      if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){      if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || ntr != sntr ){
95        newentry = true;        newentry = true;
96        OBT = L2->GetOrbitalInfo()->OBT;        OBT = L2->GetOrbitalInfo()->OBT;
97        PKT = L2->GetOrbitalInfo()->pkt_num;        PKT = L2->GetOrbitalInfo()->pkt_num;
98        atime = L2->GetOrbitalInfo()->absTime;        atime = L2->GetOrbitalInfo()->absTime;
99          sntr = ntr;
100      };      };
101    } else {    } else {
102      newentry = true;      newentry = true;
# Line 86  void CaloNuclei::Process(){ Line 104  void CaloNuclei::Process(){
104    //    //
105    if ( !newentry ) return;    if ( !newentry ) return;
106    //    //
107    //  printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);    tr = ntr;
108      //
109      if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
110    //    //
111    Clear();    Clear();
112    //    //
113    PamTrack *track = 0;    if ( debug ) printf(" Always calculate stdedx1 \n");
114    track = L2->GetTrack(0);    //
115      // Always calculate stdedx1
116    //    //
117    Int_t view = 0;    Int_t view = 0;
118    Int_t plane = 0;    Int_t plane = 0;
119    Int_t strip = 0;    Int_t strip = 0;
120      Int_t indx = 0;
121      Float_t vfpl[96];
122      Int_t stfpl[96];
123      memset(vfpl, 0, 96*sizeof(Float_t));
124      memset(stfpl, 0, 96*sizeof(Int_t));
125    Float_t mip = 0.;    Float_t mip = 0.;
126      for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
127        //
128        mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
129        //
130        // put in vfpl vector the energy release on the first plane
131        //
132        if ( strip != -1 && view == 1 && plane == 0 ) {
133          stfpl[indx] = strip;
134          vfpl[indx] = mip;
135          indx++;
136        };
137        //
138      };
139      //
140      if ( debug ) printf(" find energy released along the strip of maximum on the first plane and on the two neighbour strips \n");
141      //
142      // find energy released along the strip of maximum on the first plane and on the two neighbour strips
143      //
144      if ( indx > 0 ){
145        Int_t mindx = (Int_t)TMath::LocMax(indx,stfpl);
146        for (Int_t ii=0; ii<indx; ii++){
147          if ( stfpl[ii] == stfpl[mindx] ) stdedx1 += vfpl[ii];
148          if ( (mindx-1)>=0 && stfpl[ii] == stfpl[mindx-1] ) stdedx1 += vfpl[ii];
149          if ( (mindx+1)<96 && stfpl[ii] == stfpl[mindx+1] ) stdedx1 += vfpl[ii];
150        };
151      } else {
152        stdedx1 = 0.;
153      };
154      //
155      if ( debug ) printf(" if ( !usetrack ) return: usetrack %i ntr %i \n",usetrack,ntr);
156      //
157      //
158      //  if ( !usetrack ) return;
159      //
160      PamTrack *ptrack = 0;
161      CaloTrkVar *track = 0;
162      //
163      if ( usetrack ){
164        if ( ntr >= 0 ){
165          ptrack = L2->GetTrack(ntr);
166          if ( ptrack ) track = ptrack->GetCaloTrack();
167        } else {
168          track = L2->GetCaloStoredTrack(ntr);
169        };
170        //
171        if ( !track && ntr >= 0 ){
172          printf(" ERROR: cannot find any track!\n");
173          printf(" ERROR: CaloNuclei variables not completely filled \n");
174          return;  
175        };
176      } else {
177        if ( ntr >= 0 ){
178          if ( debug ) printf(" ERROR: you asked not to use a track but you are looking for track number %i !\n",ntr);
179          if ( debug ) printf(" ERROR: CaloNuclei variables not completely filled \n");
180          return;      
181        };
182      };
183      //
184    //  Float_t defethr = 6. * 0.90;    //  Float_t defethr = 6. * 0.90;
185    Float_t defethr = 6.25; // = (sqrt(9) - 0.5) ** 2.;    Float_t defethr = 6.25; // = (sqrt(9) - 0.5) ** 2.;
186    //    //
# Line 106  void CaloNuclei::Process(){ Line 190  void CaloNuclei::Process(){
190      //      //
191      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
192      //      //
193      if ( strip != -1 &&      if ( ntr >= 0 ){
194           view == 1   &&        //
195           plane == 0  &&        if ( strip != -1 &&
196           ( strip == (track->GetCaloTrack()->tibar[0][1]-1) || strip == (track->GetCaloTrack()->tibar[0][1]-2) || strip == (track->GetCaloTrack()->tibar[0][1]) )             view == 1   &&
197           && true ){                   plane == 0  &&
198        dedx1 += mip;             ( strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) || strip == (track->tibar[0][1]) )
199      };             && true ){      
200      if ( strip != -1 &&          dedx1 += mip;
201           (( view == 1 && ( plane == 0 || plane == 1 ) ) ||        };
202           ( view == 0 && plane == 0 )) &&        if ( strip != -1 &&
203           (( view == 0 && ( strip == track->GetCaloTrack()->tibar[0][0] || strip == (track->GetCaloTrack()->tibar[0][0]-1) || strip == (track->GetCaloTrack()->tibar[0][0]-2) )) ||             (( view == 1 && ( plane == 0 || plane == 1 ) ) ||
204            ( view == 1 && ( strip == track->GetCaloTrack()->tibar[0][1] || strip == (track->GetCaloTrack()->tibar[0][1]-1) || strip == (track->GetCaloTrack()->tibar[0][1]-2) )) ||              ( view == 0 && plane == 0 )) &&
205            ( view == 1 && ( strip == track->GetCaloTrack()->tibar[1][1] || strip == (track->GetCaloTrack()->tibar[1][1]-1) || strip == (track->GetCaloTrack()->tibar[1][1]-2) ))) &&             (( view == 0 && ( strip == track->tibar[0][0] || strip == (track->tibar[0][0]-1) || strip == (track->tibar[0][0]-2) )) ||
206           true ){              ( view == 1 && ( strip == track->tibar[0][1] || strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) )) ||
207        dedx3 += mip;              ( view == 1 && ( strip == track->tibar[1][1] || strip == (track->tibar[1][1]-1) || strip == (track->tibar[1][1]-2) ))) &&
208               true ){
209            dedx3 += mip;
210          };
211        } else {
212          //
213          if ( strip != -1 &&
214               view == 1   &&
215               plane == 0  &&
216               ( strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) || strip == (L2->GetCaloLevel2()->cibar[0][1]) )
217               && true ){      
218            dedx1 += mip;
219          };
220          if ( strip != -1 &&
221               (( view == 1 && ( plane == 0 || plane == 1 ) ) ||
222                ( view == 0 && plane == 0 )) &&
223               (( view == 0 && ( strip == L2->GetCaloLevel2()->cibar[0][0] || strip == (L2->GetCaloLevel2()->cibar[0][0]-1) || strip == (L2->GetCaloLevel2()->cibar[0][0]-2) )) ||
224                ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[0][1] || strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) )) ||
225                ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[1][1] || strip == (L2->GetCaloLevel2()->cibar[1][1]-1) || strip == (L2->GetCaloLevel2()->cibar[1][1]-2) ))) &&
226               true ){
227            dedx3 += mip;
228          };
229      };      };
230      //          //    
231    };    };
# Line 133  void CaloNuclei::Process(){ Line 238  void CaloNuclei::Process(){
238    //    //
239   retry:   retry:
240    //    //
241    //  printf("retry\n");    if ( debug ) printf("retry\n");
242    //    //
243    interplane = 0;    interplane = 0;
244    //    //
# Line 154  void CaloNuclei::Process(){ Line 259  void CaloNuclei::Process(){
259      //      //
260      mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);          mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);    
261      //      //
262      if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&      if ( ntr >= 0 ){
263           ( strip == (track->GetCaloTrack()->tibar[plane][view]-1) || strip == (track->GetCaloTrack()->tibar[plane][view]-2) || strip == (track->GetCaloTrack()->tibar[plane][view]) )        if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
264           && true ){                   ( strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) )
265        //      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->GetCaloTrack()->tibar[plane][view]-1,nhit[view],splane[view],sview[view]);             && true ){      
266        interpl[view] = plane;          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]);
267        interv[view] = view;          interpl[view] = plane;
268        if ( splane[view] != plane || sview[view] != view ){          interv[view] = view;
269          if ( nhit[view] > 1 ){          if ( splane[view] != plane || sview[view] != view ){
270            wmulthit[view] = true;            if ( nhit[view] > 1 ){
271            interpl[view] = splane[view];              wmulthit[view] = true;
272            interv[view] = sview[view];              //    if ( splane[view] == -1 ) splane[view] = 0; //
273            //      ii = L2->GetCaloLevel1()->istrip;              //    if ( sview[view] == -1 ) sview[view] = view; //
274                interpl[view] = splane[view];
275                interv[view] = sview[view];  
276          };          };
277          if ( plane > splane[view]+gapth ){            if ( plane > splane[view]+gapth ){
278            wgap[view] = true;              wgap[view] = true;
279            interpl[view] = splane[view];              //    if ( splane[view] == -1 ) splane[view] = 0;//
280            interv[view] = sview[view];              //    if ( sview[view] == -1 ) sview[view] = view; //
281            //      ii = L2->GetCaloLevel1()->istrip;              interpl[view] = splane[view];
282                interv[view] = sview[view];
283              };
284              splane[view] = plane;
285              sview[view] = view;
286              nhit[view] = 1;
287            } else {
288              nhit[view]++;
289            };
290          };
291        } else {
292          if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
293               ( strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) )
294               && true ){      
295            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]);
296            interpl[view] = plane;
297            interv[view] = view;
298            if ( splane[view] != plane || sview[view] != view ){
299              if ( nhit[view] > 1 ){
300                wmulthit[view] = true;
301                //    if ( splane[view] == -1 ) splane[view] = 0; //
302                //    if ( sview[view] == -1 ) sview[view] = view; //
303                interpl[view] = splane[view];
304                interv[view] = sview[view];  
305            };
306              if ( plane > splane[view]+gapth ){
307                wgap[view] = true;
308                //    if ( splane[view] == -1 ) splane[view] = 0;//
309                //    if ( sview[view] == -1 ) sview[view] = view; //
310                interpl[view] = splane[view];
311                interv[view] = sview[view];
312              };
313              splane[view] = plane;
314              sview[view] = view;
315              nhit[view] = 1;
316            } else {
317              nhit[view]++;
318          };          };
         splane[view] = plane;  
         sview[view] = view;  
         nhit[view] = 1;  
       } else {  
         nhit[view]++;  
319        };        };
320      };      };
321      //      //
# Line 185  void CaloNuclei::Process(){ Line 323  void CaloNuclei::Process(){
323      //          //    
324    };    };
325    //    //
326    //  printf("conversion interpl %i interv %i multhit %i interplane %i \n",interpl[0],interv[0],multhit,interplane);    if (debug ) printf("conversion interpl %i interv %i multhit %i interplane %i \n",interpl[0],interv[0],multhit,interplane);
327    Int_t winterplane[2] = {-1,-1};    Int_t winterplane[2] = {-1,-1};
328    //    //
329    for ( Int_t view = 0; view < 2; view++){    for ( Int_t view = 0; view < 2; view++){
# Line 224  void CaloNuclei::Process(){ Line 362  void CaloNuclei::Process(){
362      };      };
363    };    };
364    //      //  
365    //  printf("2conversion interpl %i interv %i multhit %i interplane %i \n",interpl[1],interv[1],multhit,interplane);    if ( debug ) printf("2conversion interpl %i interv %i multhit %i interplane %i \n",interpl[1],interv[1],multhit,interplane);
366    //  printf("3conversion winterpl0 %i winterpl1 %i \n",winterplane[0],winterplane[1]);    if ( debug ) printf("3conversion winterpl0 %i winterpl1 %i \n",winterplane[0],winterplane[1]);
367    //    //
368    Int_t ipl = 0;    Int_t ipl = 0;
369    if ( interplane > 0 ){    if ( interplane > 0 ){
# Line 236  void CaloNuclei::Process(){ Line 374  void CaloNuclei::Process(){
374      Int_t ind = -1;      Int_t ind = -1;
375      Int_t qsplane = -1;      Int_t qsplane = -1;
376      Int_t qsview = -1;      Int_t qsview = -1;
377        Int_t ind2 = -1;
378        Int_t qsplane2 = -1;
379        Int_t qsview2 = -1;
380      Float_t qme[200];      Float_t qme[200];
381      memset(qme,0,200*sizeof(Float_t));      memset(qme,0,200*sizeof(Float_t));
382        Float_t qme2[2112];
383        memset(qme2,0,2112*sizeof(Float_t));
384      //      //
385      while ( ii<L2->GetCaloLevel1()->istrip ){      while ( ii<L2->GetCaloLevel1()->istrip ){
386        //        //
# Line 253  void CaloNuclei::Process(){ Line 396  void CaloNuclei::Process(){
396            postq += mip;            postq += mip;
397          } else {          } else {
398            preq += mip;            preq += mip;
399            if (  strip == (track->GetCaloTrack()->tibar[plane][view]-1) || strip == (track->GetCaloTrack()->tibar[plane][view]-2) || strip == (track->GetCaloTrack()->tibar[plane][view]) ){            if ( ntr >= 0 ){
400              if ( qsplane != plane || qsview != view ){              if (  strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) ){
401                qsplane = plane;                if ( qsplane != plane || qsview != view ){
402                qsview = view;                  qsplane = plane;
403                ind++;                  qsview = view;
404                if ( ind > 199 ) printf(" AAAGH!! \n");                  ind++;
405                qme[ind] = 0.;                  if ( debug && ind > 199 ) printf(" AAAGH!! \n");
406                    qme[ind] = 0.;
407                  };
408                  qme[ind] += mip;
409              };              };
410              qme[ind] += mip;              for ( Int_t ns = 0; ns < R ; ns++){
411                  Int_t ms =  track->tibar[plane][view] - 1 - ns + (R - 1)/2;
412                  if ( strip == ms ){
413                    if ( qsplane2 != plane || qsview2 != view ){
414                      qsplane2 = plane;
415                      qsview2 = view;
416                      ind2++;
417                      if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
418                      qme2[ind2] = 0.;
419                    };
420                    qme2[ind2] += mip;
421                  };
422                };  
423              } else {
424                if (  strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) ){
425                  if ( qsplane != plane || qsview != view ){
426                    qsplane = plane;
427                    qsview = view;
428                    ind++;
429                    if ( debug && ind > 199 ) printf(" AAAGH!! \n");
430                    qme[ind] = 0.;
431                  };
432                  qme[ind] += mip;
433                };
434                for ( Int_t ns = 0; ns < R ; ns++){
435                  Int_t ms =  L2->GetCaloLevel2()->cibar[plane][view] - 1 - ns + (R - 1)/2;
436                  if ( strip == ms ){
437                    if ( qsplane2 != plane || qsview2 != view ){
438                      qsplane2 = plane;
439                      qsview2 = view;
440                      ind2++;
441                      if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
442                      qme2[ind2] = 0.;
443                    };
444                    qme2[ind2] += mip;
445                  };
446                };  
447            };            };
448          };                };      
449          //          //
# Line 273  void CaloNuclei::Process(){ Line 455  void CaloNuclei::Process(){
455      //      //
456      // here we must calculate qpremean, order vector qme, select 3 lowest measurements and caculate the mean...      // here we must calculate qpremean, order vector qme, select 3 lowest measurements and caculate the mean...
457      //          //    
458      //    for (Int_t l=0; l < interplane; l++){      if ( debug ){
459      //       printf(" qme[%i] = %f \n",l,qme[l]);          for (Int_t l=0; l < interplane; l++){
460      //       //      if ( qme[ind] < 6.25 ) qme[ind] = 10000.;                printf(" qme[%i] = %f \n",l,qme[l]);  
461      //     };        };
462        };
463        //
464      Long64_t work[200];      Long64_t work[200];
465      ind = 0;      ind = 0;
466      Int_t l = 0;      Int_t l = 0;
467        Int_t RN = 0;
468      Float_t qm = 0.;      Float_t qm = 0.;
469        Float_t qm2 = 0.;
470        //
471        Float_t qmt = ethr*0.8; // *0.9
472        //
473      Int_t uplim = TMath::Max(3,N);      Int_t uplim = TMath::Max(3,N);
474        //
475      while ( l < uplim && ind < interplane ){      while ( l < uplim && ind < interplane ){
476        qm = TMath::KOrdStat(interplane,qme,ind,work);        qm = TMath::KOrdStat(interplane,qme,ind,work);
477        if ( qm >= mesethr*0.9 ){        if ( qm >= qmt ){
478          if ( l < 3 ) qpremean += qm;          if ( l < 3 ){
479          if ( l < N ) qpremeanN += qm;            qpremean += qm;
480              RN++;
481            };
482          l++;          l++;
483          //      printf(" value no %i qm %f \n",l,qm);          if ( debug ) printf(" value no %i qm %f qmt %f \n",l,qm,qmt);
484        };        };
485        ind++;        ind++;
486      };      };
487      //      //
488      qpremean /= l;      qpremean /= (Float_t)RN;
     qpremeanN /= N;  
489      //      //
490      //  printf(" charge is %f \n",sqrt(qpremean));      ind = 0;
491        l = 0;
492        RN = 0;
493        while ( l < uplim && ind < interplane ){
494          qm2 = TMath::KOrdStat(interplane,qme2,ind,work);
495          if ( qm2 >= qmt ){        
496            if ( l < N ){
497              qpremeanN += qm2;
498              RN++;
499            };
500            l++;
501            if ( debug ) printf(" qm2 value no %i qm %f qmt %f RN %i \n",l,qm2,qmt,RN);
502          };
503          ind++;
504        };
505      //      //
506      //    printf("preq postq\n");          qpremeanN /= (Float_t)RN;
507        UN = RN;
508        //
509        if ( debug ) printf(" charge is %f \n",sqrt(qpremean));
510      //      //
511      if ( mesethr != ethr && interplane >= 3 && !aldone ){      if ( mesethr != ethr && interplane >= 3 && !aldone ){
512        Float_t mesethr2 = (sqrt(qpremean) - 0.50)*(sqrt(qpremean) - 0.50);        Float_t mesethr2 = (sqrt(qpremean) - 0.50)*(sqrt(qpremean) - 0.50);
# Line 308  void CaloNuclei::Process(){ Line 516  void CaloNuclei::Process(){
516          mesethr = mesethr2;          mesethr = mesethr2;
517        };        };
518        aldone = true;        aldone = true;
519        if ( mesethr > defethr ) goto retry;        if ( mesethr > defethr ){
520            interplane = 0;
521            preq = 0.;
522            postq = 0.;
523            qpremean = 0.;
524            qpremeanN = 0.;
525            multhit = false;
526            gap = false;
527            goto retry;
528          };
529      };      };
530    };    };
531    //    //
532    //  printf(" esci \n");    if ( debug ) this->Print();
533      if ( debug ) printf(" esci \n");
534    //    //
535  };  };
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.23