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

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

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

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