/[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 by mocchiut, Tue Apr 3 15:55:16 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 20  CaloNuclei::CaloNuclei(PamLevel2 *l2p){ Line 23  CaloNuclei::CaloNuclei(PamLevel2 *l2p){
23    PKT = 0;    PKT = 0;
24    atime = 0;    atime = 0;
25    N = 5;    N = 5;
26      R = 3;
27      //
28      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;
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 40  void CaloNuclei::Clear(){ Line 61  void CaloNuclei::Clear(){
61    
62  void CaloNuclei::Print(){  void CaloNuclei::Print(){
63    //    //
64    printf("==============================\n");    Process();
65    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);    //
66    printf(" interplane: %i\n",interplane);    printf("========================================================================\n");
67    printf(" ethr: %f \n",ethr);    printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack);
68    printf(" dedx1: %f \n",dedx1);    printf(" interplane [number of available dE/dx before interaction]:....... %i\n",interplane);
69    printf(" dedx3: %f \n",dedx3);    printf(" ethr [threshold used to determine interplane]:................... %f \n",ethr);
70    printf(" multhit: %i \n",multhit);    printf(" dedx1 [dE/dx from the first calorimeter plane]:.................. %f \n",dedx1);
71    printf(" gap: %i \n",gap);    printf(" stdedx1 [dE/dx from the first calorimeter plane standalone]:..... %f \n",stdedx1);
72    printf(" preq: %f \n",preq);    printf(" dedx3 [dE/dx (average) if the first 3 Si planes]:................ %f \n",dedx3);
73    printf(" postq: %f \n",postq);    printf(" multhit [true if interplane determined by multiple hits]:........ %i \n",multhit);
74    printf(" qpremean: %f \n",qpremean);    printf(" gap [true if interplane determined by a gap]:.................... %i \n",gap);
75    printf(" N: %i \n",N);    printf(" preq [total energy in MIP before the interaction plane]:......... %f \n",preq);
76    printf(" qpremeanN: %f \n",qpremeanN);    printf(" postq [total energy in MIP after the interaction plane]:......... %f \n",postq);
77    printf("==============================\n");    printf(" qpremean [truncated mean using 3 planes and 3 strips]:........... %f \n",qpremean);
78      printf(" N [no of used plane]:............................................ %i \n",N);
79      printf(" R [no strip used per plane ]:.................................... %i \n",R);
80      printf(" qpremeanN [truncated mean using N planes and R strips]:.......... %f \n",qpremeanN);
81      printf(" qNmin1 [truncated mean using N-1 planes and R strips]: .......... %f \n",qNmin1);
82      printf(" maxrel [dE/dx of strip with maximum release (I plane)]:.......... %f \n",maxrel);
83      printf(" ZCalo_maxrel_b [Z from maximum release in I Calo plane vs beta].. %f \n",ZCalo_maxrel_b);
84      printf(" ZCalo_dedx_b [Z from dedx in I Calo plane vs beta].. ............ %f \n",ZCalo_dedx_b);
85      printf(" ZCalo_dedx_defl [Z from dedx in I Calo plane vs deflection....... %f \n",ZCalo_dedx_defl);
86      printf(" ZCalo_Nmin1_defl [Z from truncated mean (N-1 pl) vs deflection].. %f \n",ZCalo_Nmin1_defl);
87      printf("========================================================================\n");
88    //    //
89  };  };
90    
# Line 64  void CaloNuclei::Delete(){ Line 95  void CaloNuclei::Delete(){
95    
96    
97  void CaloNuclei::Process(){  void CaloNuclei::Process(){
98    //    Process(0);
99    };
100    
101    void CaloNuclei::Process(Int_t ntr){
102      //  
103    if ( !L2 ){    if ( !L2 ){
104      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");
105      printf(" ERROR: CaloNuclei variables not filled \n");      printf(" ERROR: CaloNuclei variables not filled \n");
# Line 74  void CaloNuclei::Process(){ Line 109  void CaloNuclei::Process(){
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 86  void CaloNuclei::Process(){ Line 122  void CaloNuclei::Process(){
122    //    //
123    if ( !newentry ) return;    if ( !newentry ) return;
124    //    //
125    //  printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);    tr = ntr;
126      //
127      if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
128    //    //
129    Clear();    Clear();
130    //    //
131    PamTrack *track = 0;    if ( debug ) printf(" Always calculate stdedx1 \n");
   track = L2->GetTrack(0);  
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 106  void CaloNuclei::Process(){ Line 214  void CaloNuclei::Process(){
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 133  void CaloNuclei::Process(){ Line 262  void CaloNuclei::Process(){
262    //    //
263   retry:   retry:
264    //    //
265    //  printf("retry\n");    if ( debug ) printf("retry\n");
266    //    //
267    interplane = 0;    interplane = 0;
268    //    //
# Line 154  void CaloNuclei::Process(){ Line 283  void CaloNuclei::Process(){
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        //      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            interpl[view] = splane[view];              wmulthit[view] = true;
296            interv[view] = sview[view];              //    if ( splane[view] == -1 ) splane[view] = 0; //
297            //      ii = L2->GetCaloLevel1()->istrip;              //    if ( sview[view] == -1 ) sview[view] = view; //
298                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            interpl[view] = splane[view];              //    if ( splane[view] == -1 ) splane[view] = 0;//
304            interv[view] = sview[view];              //    if ( sview[view] == -1 ) sview[view] = view; //
305            //      ii = L2->GetCaloLevel1()->istrip;              interpl[view] = splane[view];
306                interv[view] = sview[view];
307              };
308              splane[view] = plane;
309              sview[view] = view;
310              nhit[view] = 1;
311            } else {
312              nhit[view]++;
313            };
314          };
315        } else {
316          if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
317               ( strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) )
318               && true ){      
319            if ( debug ) printf(" inside loop: ii %i mip %f view %i plane %i strip %i cibar %i nhit %i splane %i sview %i \n",ii,mip,view,plane,strip,L2->GetCaloLevel2()->cibar[plane][view]-1,nhit[view],splane[view],sview[view]);
320            interpl[view] = plane;
321            interv[view] = view;
322            if ( splane[view] != plane || sview[view] != view ){
323              if ( nhit[view] > 1 ){
324                wmulthit[view] = true;
325                //    if ( splane[view] == -1 ) splane[view] = 0; //
326                //    if ( sview[view] == -1 ) sview[view] = view; //
327                interpl[view] = splane[view];
328                interv[view] = sview[view];  
329            };
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 185  void CaloNuclei::Process(){ Line 347  void CaloNuclei::Process(){
347      //          //    
348    };    };
349    //    //
350    //  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);
351    Int_t winterplane[2] = {-1,-1};    Int_t winterplane[2] = {-1,-1};
352    //    //
353    for ( Int_t view = 0; view < 2; view++){    for ( Int_t view = 0; view < 2; view++){
# Line 224  void CaloNuclei::Process(){ Line 386  void CaloNuclei::Process(){
386      };      };
387    };    };
388    //      //  
389    //  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);
390    //  printf("3conversion winterpl0 %i winterpl1 %i \n",winterplane[0],winterplane[1]);    if ( debug ) printf("3conversion winterpl0 %i winterpl1 %i \n",winterplane[0],winterplane[1]);
391    //    //
392    Int_t ipl = 0;    Int_t ipl = 0;
393    if ( interplane > 0 ){    if ( interplane > 0 ){
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;
401      Int_t qsview = -1;      Int_t qsview = -1;
402        Int_t ind2 = -1;
403        Int_t qsplane2 = -1;
404        Int_t qsview2 = -1;
405      Float_t qme[200];      Float_t qme[200];
406      memset(qme,0,200*sizeof(Float_t));      memset(qme,0,200*sizeof(Float_t));
407        Float_t qme2[2112];
408        memset(qme2,0,2112*sizeof(Float_t));
409      //      //
410      while ( ii<L2->GetCaloLevel1()->istrip ){      while ( ii<L2->GetCaloLevel1()->istrip ){
411        //        //
# Line 253  void CaloNuclei::Process(){ Line 421  void CaloNuclei::Process(){
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 ( 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                  if ( strip == ms ){
438                    if ( qsplane2 != plane || qsview2 != view ){
439                      qsplane2 = plane;
440                      qsview2 = view;
441                      ind2++;
442                      if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
443                      qme2[ind2] = 0.;
444                    };
445                    qme2[ind2] += mip;
446                  };
447                };  
448              } else {
449                if (  strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) ){
450                  if ( qsplane != plane || qsview != view ){
451                    qsplane = plane;
452                    qsview = view;
453                    ind++;
454                    if ( debug && ind > 199 ) printf(" AAAGH!! \n");
455                    qme[ind] = 0.;
456                  };
457                  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          //          //
# Line 270  void CaloNuclei::Process(){ Line 477  void CaloNuclei::Process(){
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      //    for (Int_t l=0; l < interplane; l++){      if ( debug ){
486      //       printf(" qme[%i] = %f \n",l,qme[l]);          for (Int_t l=0; l < interplane; l++){
487      //       //      if ( qme[ind] < 6.25 ) qme[ind] = 10000.;                printf(" qme[%i] = %f \n",l,qme[l]);  
488      //     };        };
489        };
490        //
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.;
497        //
498        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 >= mesethr*0.9 ){        if ( qm >= qmt ){
506          if ( l < 3 ) qpremean += qm;          if ( l < 3 ){
507          if ( l < N ) qpremeanN += qm;            qpremean += qm;
508              RN++;
509            };
510          l++;          l++;
511          //      printf(" value no %i qm %f \n",l,qm);          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      qpremeanN /= N;      ind = 0;
518      //      l = 0;
519      //  printf(" charge is %f \n",sqrt(qpremean));      RN = 0;
520        while ( l < uplim && ind < interplane ){
521          qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
522          if ( qm2 >= qmt ){        
523            if ( l < N ){
524              qpremeanN += qm2;
525              RN++;
526            };
527            l++;
528            if ( debug ) printf(" qm2 value no %i qm %f qmt %f RN %i \n",l,qm2,qmt,RN);
529          };
530          ind++;
531        };
532        ////////////////////////////////////
533        //to calculate qNmin1///////////////
534        ///////////////////////////////////
535        //values under threshold
536        qm2=0;
537        ind = 0;
538        l = 0;
539        RN = 0;
540        S2=0;
541        while ( l < uplim2 && ind<interplane){
542          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      //    printf("preq postq\n");          if ( debug ) printf(" charge is %f \n",sqrt(qpremean));
577      //      //
578      if ( mesethr != ethr && interplane >= 3 && !aldone ){      if ( mesethr != ethr && interplane >= 3 && !aldone ){
579        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 583  void CaloNuclei::Process(){
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    //  printf(" esci \n");    if ( debug ) this->Print();
872      if ( debug ) printf(" esci \n");
873    //    //
874  };  };
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.23