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