/[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.16 by mocchiut, Fri Nov 28 16:01:26 2008 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    //    //
32  };  };
33    
34  void CaloNuclei::Clear(){  void CaloNuclei::Clear(){
35    //    //
36      UN = 0;
37    tr = 0;    tr = 0;
38      sntr = 0;
39    interplane = 0;    interplane = 0;
40    preq = 0.;    preq = 0.;
41    postq = 0.;    postq = 0.;
42      stdedx1 = 0.;
43      ethr = 0.;
44    dedx1 = 0.;    dedx1 = 0.;
45    dedx3 = 0.;    dedx3 = 0.;
46    qpremean = 0.;    qpremean = 0.;
47    qpremeanN = 0.;    qpremeanN = 0.;
48      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    //    //
57    multhit = false;    multhit = false;
58    gap = false;    gap = false;
# Line 47  void CaloNuclei::Print(){ Line 64  void CaloNuclei::Print(){
64    Process();    Process();
65    //    //
66    printf("========================================================================\n");    printf("========================================================================\n");
67    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);
68    printf(" interplane [number of available dE/dx before interaction]: %i\n",interplane);    printf(" interplane [number of available dE/dx before interaction]:....... %i\n",interplane);
69    printf(" ethr [threshold used to determine interplane]:............ %f \n",ethr);    printf(" ethr [threshold used to determine interplane]:................... %f \n",ethr);
70    printf(" dedx1 [dE/dx from the first calorimeter plane]:........... %f \n",dedx1);    printf(" dedx1 [dE/dx from the first calorimeter plane]:.................. %f \n",dedx1);
71    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);
72    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);
73    printf(" gap [true if interplane determined by a gap]:............. %i \n",gap);    printf(" multhit [true if interplane determined by multiple hits]:........ %i \n",multhit);
74    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);
75    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);
76    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);
77    printf(" N [no of used plane]:..................................... %i \n",N);    printf(" qpremean [truncated mean using 3 planes and 3 strips]:........... %f \n",qpremean);
78    printf(" R [no strip used per plane ]:............................. %i \n",R);    printf(" N [no of used plane]:............................................ %i \n",N);
79    printf(" qpremeanN [truncated mean using N planes and R strips]:... %f \n",qpremeanN);    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    printf("========================================================================\n");    printf("========================================================================\n");
88    //    //
89  };  };
# Line 85  void CaloNuclei::Process(Int_t ntr){ Line 109  void CaloNuclei::Process(Int_t ntr){
109    Bool_t newentry = false;    Bool_t newentry = false;
110    //    //
111    if ( L2->IsORB() ){    if ( L2->IsORB() ){
112      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 ){
113        newentry = true;        newentry = true;
114        OBT = L2->GetOrbitalInfo()->OBT;        OBT = L2->GetOrbitalInfo()->OBT;
115        PKT = L2->GetOrbitalInfo()->pkt_num;        PKT = L2->GetOrbitalInfo()->pkt_num;
116        atime = L2->GetOrbitalInfo()->absTime;        atime = L2->GetOrbitalInfo()->absTime;
117          sntr = ntr;
118      };      };
119    } else {    } else {
120      newentry = true;      newentry = true;
# Line 103  void CaloNuclei::Process(Int_t ntr){ Line 128  void CaloNuclei::Process(Int_t ntr){
128    //    //
129    Clear();    Clear();
130    //    //
131    PamTrack *track = 0;    if ( debug ) printf(" Always calculate stdedx1 \n");
   track = L2->GetTrack(ntr);  
132    //    //
133      // Always calculate stdedx1 and maxrel
134      //
135      Int_t cont=0;
136    Int_t view = 0;    Int_t view = 0;
137    Int_t plane = 0;    Int_t plane = 0;
138    Int_t strip = 0;    Int_t strip = 0;
139      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    Float_t mip = 0.;    Float_t mip = 0.;
145      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      if ( debug ) printf(" find energy released along the strip of maximum on the first plane and on the two neighbour strips \n");
160      //
161      // find energy released along the strip of maximum on the first plane and on the two neighbour strips
162      //
163      if ( indx > 0 ){
164        Int_t mindx = (Int_t)TMath::LocMax(indx,vfpl);
165        for (Int_t ii=0; ii<indx; ii++){
166          if ( stfpl[ii] == stfpl[mindx] ) stdedx1 += vfpl[ii];
167          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        };
172      maxrel = vfpl[mindx];
173      } else {
174        stdedx1 = 0.;
175        maxrel = 0.;
176      };
177      // cout<<stdedx1<<"      "<<maxrel<<"\n";
178      //
179      if ( debug ) printf(" if ( !usetrack ) return: usetrack %i ntr %i \n",usetrack,ntr);
180      //
181      //
182      //  if ( !usetrack ) return;
183      //
184      PamTrack *ptrack = 0;
185      CaloTrkVar *track = 0;
186      //
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      } else {
201        if ( ntr >= 0 ){
202          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          return;      
205        };
206      };
207      //
208    //  Float_t defethr = 6. * 0.90;    //  Float_t defethr = 6. * 0.90;
209    Float_t defethr = 6.25; // = (sqrt(9) - 0.5) ** 2.;    Float_t defethr = 6.25; // = (sqrt(9) - 0.5) ** 2.;
210    //    //
# Line 119  void CaloNuclei::Process(Int_t ntr){ Line 214  void CaloNuclei::Process(Int_t ntr){
214      //      //
215      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
216      //      //
217      if ( strip != -1 &&      if ( ntr >= 0 ){
218           view == 1   &&        //
219           plane == 0  &&        if ( strip != -1 &&
220           ( strip == (track->GetCaloTrack()->tibar[0][1]-1) || strip == (track->GetCaloTrack()->tibar[0][1]-2) || strip == (track->GetCaloTrack()->tibar[0][1]) )             view == 1   &&
221           && true ){                   plane == 0  &&
222        dedx1 += mip;             ( strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) || strip == (track->tibar[0][1]) )
223      };             && true ){      
224      if ( strip != -1 &&          dedx1 += mip;
225           (( view == 1 && ( plane == 0 || plane == 1 ) ) ||        };
226           ( view == 0 && plane == 0 )) &&        if ( strip != -1 &&
227           (( 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 ) ) ||
228            ( 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 )) &&
229            ( 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) )) ||
230           true ){              ( view == 1 && ( strip == track->tibar[0][1] || strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) )) ||
231        dedx3 += mip;              ( 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      };      };
254      //          //    
255    };    };
# Line 167  void CaloNuclei::Process(Int_t ntr){ Line 283  void CaloNuclei::Process(Int_t ntr){
283      //      //
284      mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);          mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);    
285      //      //
286      if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&      if ( ntr >= 0 ){
287           ( 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] &&
288           && true ){                   ( strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) )
289        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]);             && true ){      
290        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]);
291        interv[view] = view;          interpl[view] = plane;
292        if ( splane[view] != plane || sview[view] != view ){          interv[view] = view;
293          if ( nhit[view] > 1 ){          if ( splane[view] != plane || sview[view] != view ){
294            wmulthit[view] = true;            if ( nhit[view] > 1 ){
295            //      if ( splane[view] == -1 ) splane[view] = 0; //              wmulthit[view] = true;
296            //      if ( sview[view] == -1 ) sview[view] = view; //              //    if ( splane[view] == -1 ) splane[view] = 0; //
297            interpl[view] = splane[view];              //    if ( sview[view] == -1 ) sview[view] = view; //
298            interv[view] = sview[view];                  interpl[view] = splane[view];
299                interv[view] = sview[view];  
300          };          };
301          if ( plane > splane[view]+gapth ){            if ( plane > splane[view]+gapth ){
302            wgap[view] = true;              wgap[view] = true;
303            //      if ( splane[view] == -1 ) splane[view] = 0;//              //    if ( splane[view] == -1 ) splane[view] = 0;//
304            //      if ( sview[view] == -1 ) sview[view] = view; //              //    if ( sview[view] == -1 ) sview[view] = view; //
305            interpl[view] = splane[view];              interpl[view] = splane[view];
306            interv[view] = sview[view];              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            };
330              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          };          };
         splane[view] = plane;  
         sview[view] = view;  
         nhit[view] = 1;  
       } else {  
         nhit[view]++;  
343        };        };
344      };      };
345      //      //
# Line 247  void CaloNuclei::Process(Int_t ntr){ Line 394  void CaloNuclei::Process(Int_t ntr){
394      //      //
395      // Calculate preq, postq, qpremean      // Calculate preq, postq, qpremean
396      //      //
397        cont++;
398      ii = 0;      ii = 0;
399      Int_t ind = -1;      Int_t ind = -1;
400      Int_t qsplane = -1;      Int_t qsplane = -1;
# Line 273  void CaloNuclei::Process(Int_t ntr){ Line 421  void CaloNuclei::Process(Int_t ntr){
421            postq += mip;            postq += mip;
422          } else {          } else {
423            preq += mip;            preq += mip;
424            if (  strip == (track->GetCaloTrack()->tibar[plane][view]-1) || strip == (track->GetCaloTrack()->tibar[plane][view]-2) || strip == (track->GetCaloTrack()->tibar[plane][view]) ){            if ( ntr >= 0 ){
425              if ( qsplane != plane || qsview != view ){              if (  strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) ){
426                qsplane = plane;                if ( qsplane != plane || qsview != view ){
427                qsview = view;                  qsplane = plane;
428                ind++;                  qsview = view;
429                if ( debug && ind > 199 ) printf(" AAAGH!! \n");                  ind++;
430                qme[ind] = 0.;                  if ( debug && ind > 199 ) printf(" AAAGH!! \n");
431                    qme[ind] = 0.;
432                  };
433                  qme[ind] += mip;
434              };              };
435              qme[ind] += mip;              for ( Int_t ns = 0; ns < R ; ns++){
436            };                Int_t ms =  track->tibar[plane][view] - 1 - ns + (R - 1)/2;
437            for ( Int_t ns = 0; ns < R ; ns++){                if ( strip == ms ){
438              Int_t ms =  track->GetCaloTrack()->tibar[plane][view] - 1 - ns + (R - 1)/2;                  if ( qsplane2 != plane || qsview2 != view ){
439              if ( strip == ms ){                    qsplane2 = plane;
440                if ( qsplane2 != plane || qsview2 != view ){                    qsview2 = view;
441                  qsplane2 = plane;                    ind2++;
442                  qsview2 = view;                    if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
443                  ind2++;                    qme2[ind2] = 0.;
444                  if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");                  };
445                  qme2[ind2] = 0.;                  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                };                };
457                qme2[ind2] += mip;                qme[ind] += mip;
458              };              };
459            };                  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          };                };      
474          //          //
475        };        };
# Line 303  void CaloNuclei::Process(Int_t ntr){ Line 477  void CaloNuclei::Process(Int_t ntr){
477        ii++;        ii++;
478        //            //    
479      };      };
480      //  
481    
482     //
483      // 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...
484      //          //    
485      if ( debug ){      if ( debug ){
# Line 315  void CaloNuclei::Process(Int_t ntr){ Line 491  void CaloNuclei::Process(Int_t ntr){
491      Long64_t work[200];      Long64_t work[200];
492      ind = 0;      ind = 0;
493      Int_t l = 0;      Int_t l = 0;
494        Int_t RN = 0;
495      Float_t qm = 0.;      Float_t qm = 0.;
496      Float_t qm2 = 0.;      Float_t qm2 = 0.;
497      //      //
498      Float_t qmt = ethr*0.8; // *0.9      Float_t qmt = ethr*0.8; // *0.9
499      //      //
500      Int_t uplim = TMath::Max(3,N);      Int_t uplim = TMath::Max(3,N);
501        Int_t uplim2 = interplane-1;
502      //      //
503      while ( l < uplim && ind < interplane ){      while ( l < uplim && ind < interplane ){
504        qm = TMath::KOrdStat(interplane,qme,ind,work);        qm = TMath::KOrdStat((Long64_t)interplane,qme,(Long64_t)ind,work);
505        if ( qm >= qmt ){        if ( qm >= qmt ){
506          if ( l < 3 ) qpremean += qm;          if ( l < 3 ){
507              qpremean += qm;
508              RN++;
509            };
510          l++;          l++;
511          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);
512        };        };
513        ind++;        ind++;
514      };      };
515      //      //
516      qpremean /= l;      qpremean /= (Float_t)RN;
     //  
517      ind = 0;      ind = 0;
518      l = 0;      l = 0;
519        RN = 0;
520      while ( l < uplim && ind < interplane ){      while ( l < uplim && ind < interplane ){
521        qm2 = TMath::KOrdStat(interplane,qme2,ind,work);        qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
522        if ( qm2 >= qmt ){                if ( qm2 >= qmt ){        
523          if ( l < N ) qpremeanN += qm2;          if ( l < N ){
524              qpremeanN += qm2;
525              RN++;
526            };
527          l++;          l++;
528          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);
529        };        };
530        ind++;        ind++;
531      };      };
532      //      ////////////////////////////////////
533      qpremeanN /= l;      //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          qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
543          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          qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
552          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        qpremeanN /= (Float_t)RN;
563        qNmin1_w /= (Float_t)RN;
564        UN = RN;
565        ///////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      //      //
576      if ( debug ) printf(" charge is %f \n",sqrt(qpremean));      if ( debug ) printf(" charge is %f \n",sqrt(qpremean));
577      //      //
# Line 358  void CaloNuclei::Process(Int_t ntr){ Line 583  void CaloNuclei::Process(Int_t ntr){
583          mesethr = mesethr2;          mesethr = mesethr2;
584        };        };
585        aldone = true;        aldone = true;
586        if ( mesethr > defethr ) goto retry;        if ( mesethr > defethr ){
587            interplane = 0;
588            preq = 0.;
589            postq = 0.;
590            qpremean = 0.;
591            qpremeanN = 0.;
592            qNmin1 = 0;
593            multhit = false;
594            gap = false;
595            goto retry;
596          };
597      };      };
598    };    };
599    
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    // //=======================================================================
669        
670        Float_t D0[7] = {0, 3 , 4 , 5 , 6, 8, 90};
671        Float_t E1[7] = {0 ,923.553 , 659.842, 1113.97, 3037.25, 3034.84, 0};
672        Float_t E2[7] = {0 ,6.92574 , 5.08865, 5.29349, 6.41442, 5.52969, 0};
673        Float_t E3[7] = {0 ,9.7227  , 13.18, 23.5444, 38.2057, 63.6784, 80000};
674    
675        Float_t xx1[7],yy1[7];
676        n1 = 7;
677        
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    // =======================================================================
718        
719        Float_t F0[7] = {0.,3. ,4., 5. , 6., 8, 90};
720        Float_t G1[7] = {0 ,642.935 , 848.684, 1346.05, 3238.82, 3468.6, 0};
721        Float_t G2[7] = {0 ,6.2038 , 5.51723, 5.65265, 6.54089, 5.72723, 0};
722        Float_t G3[7] = {0 ,9.2421 , 13.9858, 25.3912, 39.6332, 64.5674, 80000};
723    
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    //=======================================================================
767        
768        //new
769         Float_t H0[7] = {0, 3 , 4 , 5 , 6, 8, 90 };
770         Float_t I1[7] = {0 , 56.1019, 101.673, 109.225, 150.599, 388.531, 0};
771         Float_t I2[7] = {0 , -12.5581, -22.5543, -15.9823, -28.2207, -93.6871, 0};
772         Float_t I3[7] = {0 , 11.6218, 19.664, 32.1817, 45.7527, 84.5992, 80000};
773        
774    
775    //     Float_t H0[7] = {0, 3 , 4 , 5 , 6, 8, 90 };
776    //     Float_t I1[7] = {0 , 56.1019, 101.673, 155, 150.599, 388.531, 0};
777    //     Float_t I2[7] = {0 , -12.5581, -22.5543, -35.6217, -28.2207, -93.6871, 0};
778    //     Float_t I3[7] = {0 , 11.6218, 19.664, 34.3311, 45.7527, 84.5992, 8000};
779    
780    
781        
782        charge = 1000.;
783        mip=0;
784        Float_t defl=0;
785        
786    
787        if (beta<2.) { // it makes no sense to allow beta=5 or so...
788    
789          if( L2->GetTrkLevel2()->GetNTracks()>=1 ){
790            PamTrack *TRKtrack = L2->GetTrack(0);
791            mip=dedx1;
792            if (mip==0) mip=stdedx1;
793            defl=TRKtrack->GetTrkTrack()->al[4];
794            
795            
796            if (mip>0 && defl<0.7 && defl>0)   {
797              
798              Float_t ym = mip;
799              Float_t xb = defl;
800              
801              for ( Int_t jj=0; jj<n1; jj++ ){
802                xx1[jj] = I1[jj]*xb*xb+I2[jj]*xb+I3[jj];
803                yy1[jj] = H0[jj]*H0[jj] ;
804              }
805              
806              TGraph *gr1 = new TGraph(n1,xx1,yy1);
807              TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
808              Float_t chelp = spl1->Eval(ym);
809              charge = TMath::Sqrt(chelp);
810              gr1->Delete();
811              spl1->Delete();
812              
813            } // if (mip1>0 && defl<0.5 && defl>0)
814          }//Ntrack>=1
815        } //beta < 100
816    
817        ZCalo_dedx_defl = charge;
818    
819    //=======================  end charge Rome: dedx vs defl  ===========================
820    
821    
822    //============================================================================================
823    //===========    charge determination Truncated mean (N-1 planes) vs. defl ===================
824    //================================     Rome method    ========================================
825    //============================================================================================
826        
827        Float_t L0[7] = {0, 3 , 4 , 5 , 6, 8, 90};
828        Float_t M1[7] = {0 , 63.0145, 120.504, 173.663, 245.33, 236.517, 0};
829        Float_t M2[7] = {0 , -15.005, -31.0635, -39.4988, -60.5011, -46.3992, 0};
830        Float_t M3[7] = {0 , 12.5037, 22.8652, 35.2907, 51.4678, 86.4155, 8000};
831    
832        charge = 1000.;
833        mip=0;
834      
835    
836        if (beta<2.) { // it makes no sense to allow beta=5 or so...
837    
838          if( L2->GetTrkLevel2()->GetNTracks()>=1 ){
839            mip=qNmin1;
840            
841            if (mip>0 && defl<0.7 && defl>0)   {
842              
843              Float_t ym = mip;
844              Float_t xb = defl;
845              
846              for ( Int_t jj=0; jj<n1; jj++ ){
847                xx1[jj] = M1[jj]*xb*xb+M2[jj]*xb+M3[jj];
848                yy1[jj] = L0[jj]*L0[jj] ;
849              }
850              
851              TGraph *gr1 = new TGraph(n1,xx1,yy1);
852              TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
853              Float_t chelp = spl1->Eval(ym);
854              charge = TMath::Sqrt(chelp);
855              gr1->Delete();
856              spl1->Delete();
857              
858            } // if (mip1>0 && defl<0.5 && defl>0)
859          }//Ntrack>=1
860        } //beta < 100
861        
862        ZCalo_Nmin1_defl = charge;
863    
864    //=======================  end charge Rome: Nmin1 vs defl  ===========================
865    
866    
867    
868    
869    
870    //    //
871      if ( debug ) this->Print();
872    if ( debug ) printf(" esci \n");    if ( debug ) printf(" esci \n");
873    //    //
874  };  };
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.23