/[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.13 by mocchiut, Fri Mar 14 09:24:35 2008 UTC
# Line 23  CaloNuclei::CaloNuclei(PamLevel2 *l2p){ Line 23  CaloNuclei::CaloNuclei(PamLevel2 *l2p){
23    R = 3;    R = 3;
24    //    //
25    debug = false;    debug = false;
26      usetrack = true;
27    //    //
28  };  };
29    
30  void CaloNuclei::Clear(){  void CaloNuclei::Clear(){
31    //    //
32      UN = 0;
33    tr = 0;    tr = 0;
34      sntr = 0;
35    interplane = 0;    interplane = 0;
36    preq = 0.;    preq = 0.;
37    postq = 0.;    postq = 0.;
38      stdedx1 = 0.;
39      ethr = 0.;
40    dedx1 = 0.;    dedx1 = 0.;
41    dedx3 = 0.;    dedx3 = 0.;
42    qpremean = 0.;    qpremean = 0.;
# Line 47  void CaloNuclei::Print(){ Line 52  void CaloNuclei::Print(){
52    Process();    Process();
53    //    //
54    printf("========================================================================\n");    printf("========================================================================\n");
55    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);
56    printf(" interplane [number of available dE/dx before interaction]: %i\n",interplane);    printf(" interplane [number of available dE/dx before interaction]:.. %i\n",interplane);
57    printf(" ethr [threshold used to determine interplane]:............ %f \n",ethr);    printf(" ethr [threshold used to determine interplane]:.............. %f \n",ethr);
58    printf(" dedx1 [dE/dx from the first calorimeter plane]:........... %f \n",dedx1);    printf(" dedx1 [dE/dx from the first calorimeter plane]:............. %f \n",dedx1);
59    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);
60    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);
61    printf(" gap [true if interplane determined by a gap]:............. %i \n",gap);    printf(" multhit [true if interplane determined by multiple hits]:... %i \n",multhit);
62    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);
63    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);
64    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);
65    printf(" N [no of used plane]:..................................... %i \n",N);    printf(" qpremean [truncated mean using 3 planes and 3 strips]:...... %f \n",qpremean);
66    printf(" R [no strip used per plane ]:............................. %i \n",R);    printf(" N [no of used plane]:....................................... %i \n",N);
67    printf(" qpremeanN [truncated mean using N planes and R strips]:... %f \n",qpremeanN);    printf(" R [no strip used per plane ]:............................... %i \n",R);
68      printf(" qpremeanN [truncated mean using N planes and R strips]:..... %f \n",qpremeanN);
69    printf("========================================================================\n");    printf("========================================================================\n");
70    //    //
71  };  };
# Line 85  void CaloNuclei::Process(Int_t ntr){ Line 91  void CaloNuclei::Process(Int_t ntr){
91    Bool_t newentry = false;    Bool_t newentry = false;
92    //    //
93    if ( L2->IsORB() ){    if ( L2->IsORB() ){
94      if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){      if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || ntr != sntr ){
95        newentry = true;        newentry = true;
96        OBT = L2->GetOrbitalInfo()->OBT;        OBT = L2->GetOrbitalInfo()->OBT;
97        PKT = L2->GetOrbitalInfo()->pkt_num;        PKT = L2->GetOrbitalInfo()->pkt_num;
98        atime = L2->GetOrbitalInfo()->absTime;        atime = L2->GetOrbitalInfo()->absTime;
99          sntr = ntr;
100      };      };
101    } else {    } else {
102      newentry = true;      newentry = true;
# Line 103  void CaloNuclei::Process(Int_t ntr){ Line 110  void CaloNuclei::Process(Int_t ntr){
110    //    //
111    Clear();    Clear();
112    //    //
113    PamTrack *track = 0;    if ( debug ) printf(" Always calculate stdedx1 \n");
114    track = L2->GetTrack(ntr);    //
115      // Always calculate stdedx1
116    //    //
117    Int_t view = 0;    Int_t view = 0;
118    Int_t plane = 0;    Int_t plane = 0;
119    Int_t strip = 0;    Int_t strip = 0;
120      Int_t indx = 0;
121      Float_t vfpl[96];
122      Int_t stfpl[96];
123      memset(vfpl, 0, 96*sizeof(Float_t));
124      memset(stfpl, 0, 96*sizeof(Int_t));
125    Float_t mip = 0.;    Float_t mip = 0.;
126      for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
127        //
128        mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
129        //
130        // put in vfpl vector the energy release on the first plane
131        //
132        if ( strip != -1 && view == 1 && plane == 0 ) {
133          stfpl[indx] = strip;
134          vfpl[indx] = mip;
135          indx++;
136        };
137        //
138      };
139      //
140      if ( debug ) printf(" find energy released along the strip of maximum on the first plane and on the two neighbour strips \n");
141      //
142      // find energy released along the strip of maximum on the first plane and on the two neighbour strips
143      //
144      if ( indx > 0 ){
145        Int_t mindx = (Int_t)TMath::LocMax(indx,vfpl);
146        for (Int_t ii=0; ii<indx; ii++){
147          if ( stfpl[ii] == stfpl[mindx] ) stdedx1 += vfpl[ii];
148          if ( (mindx-1)>=0 && stfpl[ii] == (stfpl[mindx]-1) ) stdedx1 += vfpl[ii];
149          if ( (mindx+1)<96 && stfpl[ii] == (stfpl[mindx]+1) ) stdedx1 += vfpl[ii];
150          //      if ( (mindx-1)>=0 && stfpl[ii] == stfpl[mindx-1] ) stdedx1 += vfpl[ii];
151          //      if ( (mindx+1)<96 && stfpl[ii] == stfpl[mindx+1] ) stdedx1 += vfpl[ii];
152        };
153      } else {
154        stdedx1 = 0.;
155      };
156      //
157      if ( debug ) printf(" if ( !usetrack ) return: usetrack %i ntr %i \n",usetrack,ntr);
158      //
159      //
160      //  if ( !usetrack ) return;
161      //
162      PamTrack *ptrack = 0;
163      CaloTrkVar *track = 0;
164      //
165      if ( usetrack ){
166        if ( ntr >= 0 ){
167          ptrack = L2->GetTrack(ntr);
168          if ( ptrack ) track = ptrack->GetCaloTrack();
169        } else {
170          track = L2->GetCaloStoredTrack(ntr);
171        };
172        //
173        if ( !track && ntr >= 0 ){
174          printf(" ERROR: cannot find any track!\n");
175          printf(" ERROR: CaloNuclei variables not completely filled \n");
176          return;  
177        };
178      } else {
179        if ( ntr >= 0 ){
180          if ( debug ) printf(" ERROR: you asked not to use a track but you are looking for track number %i !\n",ntr);
181          if ( debug ) printf(" ERROR: CaloNuclei variables not completely filled \n");
182          return;      
183        };
184      };
185      //
186    //  Float_t defethr = 6. * 0.90;    //  Float_t defethr = 6. * 0.90;
187    Float_t defethr = 6.25; // = (sqrt(9) - 0.5) ** 2.;    Float_t defethr = 6.25; // = (sqrt(9) - 0.5) ** 2.;
188    //    //
# Line 119  void CaloNuclei::Process(Int_t ntr){ Line 192  void CaloNuclei::Process(Int_t ntr){
192      //      //
193      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
194      //      //
195      if ( strip != -1 &&      if ( ntr >= 0 ){
196           view == 1   &&        //
197           plane == 0  &&        if ( strip != -1 &&
198           ( strip == (track->GetCaloTrack()->tibar[0][1]-1) || strip == (track->GetCaloTrack()->tibar[0][1]-2) || strip == (track->GetCaloTrack()->tibar[0][1]) )             view == 1   &&
199           && true ){                   plane == 0  &&
200        dedx1 += mip;             ( strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) || strip == (track->tibar[0][1]) )
201      };             && true ){      
202      if ( strip != -1 &&          dedx1 += mip;
203           (( view == 1 && ( plane == 0 || plane == 1 ) ) ||        };
204           ( view == 0 && plane == 0 )) &&        if ( strip != -1 &&
205           (( 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 ) ) ||
206            ( 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 )) &&
207            ( 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) )) ||
208           true ){              ( view == 1 && ( strip == track->tibar[0][1] || strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) )) ||
209        dedx3 += mip;              ( view == 1 && ( strip == track->tibar[1][1] || strip == (track->tibar[1][1]-1) || strip == (track->tibar[1][1]-2) ))) &&
210               true ){
211            dedx3 += mip;
212          };
213        } else {
214          //
215          if ( strip != -1 &&
216               view == 1   &&
217               plane == 0  &&
218               ( strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) || strip == (L2->GetCaloLevel2()->cibar[0][1]) )
219               && true ){      
220            dedx1 += mip;
221          };
222          if ( strip != -1 &&
223               (( view == 1 && ( plane == 0 || plane == 1 ) ) ||
224                ( view == 0 && plane == 0 )) &&
225               (( view == 0 && ( strip == L2->GetCaloLevel2()->cibar[0][0] || strip == (L2->GetCaloLevel2()->cibar[0][0]-1) || strip == (L2->GetCaloLevel2()->cibar[0][0]-2) )) ||
226                ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[0][1] || strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) )) ||
227                ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[1][1] || strip == (L2->GetCaloLevel2()->cibar[1][1]-1) || strip == (L2->GetCaloLevel2()->cibar[1][1]-2) ))) &&
228               true ){
229            dedx3 += mip;
230          };
231      };      };
232      //          //    
233    };    };
# Line 167  void CaloNuclei::Process(Int_t ntr){ Line 261  void CaloNuclei::Process(Int_t ntr){
261      //      //
262      mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);          mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);    
263      //      //
264      if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&      if ( ntr >= 0 ){
265           ( 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] &&
266           && true ){                   ( strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) )
267        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 ){      
268        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]);
269        interv[view] = view;          interpl[view] = plane;
270        if ( splane[view] != plane || sview[view] != view ){          interv[view] = view;
271          if ( nhit[view] > 1 ){          if ( splane[view] != plane || sview[view] != view ){
272            wmulthit[view] = true;            if ( nhit[view] > 1 ){
273            //      if ( splane[view] == -1 ) splane[view] = 0; //              wmulthit[view] = true;
274            //      if ( sview[view] == -1 ) sview[view] = view; //              //    if ( splane[view] == -1 ) splane[view] = 0; //
275            interpl[view] = splane[view];              //    if ( sview[view] == -1 ) sview[view] = view; //
276            interv[view] = sview[view];                  interpl[view] = splane[view];
277                interv[view] = sview[view];  
278          };          };
279          if ( plane > splane[view]+gapth ){            if ( plane > splane[view]+gapth ){
280            wgap[view] = true;              wgap[view] = true;
281            //      if ( splane[view] == -1 ) splane[view] = 0;//              //    if ( splane[view] == -1 ) splane[view] = 0;//
282            //      if ( sview[view] == -1 ) sview[view] = view; //              //    if ( sview[view] == -1 ) sview[view] = view; //
283            interpl[view] = splane[view];              interpl[view] = splane[view];
284            interv[view] = sview[view];              interv[view] = sview[view];
285              };
286              splane[view] = plane;
287              sview[view] = view;
288              nhit[view] = 1;
289            } else {
290              nhit[view]++;
291            };
292          };
293        } else {
294          if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
295               ( strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) )
296               && true ){      
297            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]);
298            interpl[view] = plane;
299            interv[view] = view;
300            if ( splane[view] != plane || sview[view] != view ){
301              if ( nhit[view] > 1 ){
302                wmulthit[view] = true;
303                //    if ( splane[view] == -1 ) splane[view] = 0; //
304                //    if ( sview[view] == -1 ) sview[view] = view; //
305                interpl[view] = splane[view];
306                interv[view] = sview[view];  
307            };
308              if ( plane > splane[view]+gapth ){
309                wgap[view] = true;
310                //    if ( splane[view] == -1 ) splane[view] = 0;//
311                //    if ( sview[view] == -1 ) sview[view] = view; //
312                interpl[view] = splane[view];
313                interv[view] = sview[view];
314              };
315              splane[view] = plane;
316              sview[view] = view;
317              nhit[view] = 1;
318            } else {
319              nhit[view]++;
320          };          };
         splane[view] = plane;  
         sview[view] = view;  
         nhit[view] = 1;  
       } else {  
         nhit[view]++;  
321        };        };
322      };      };
323      //      //
# Line 273  void CaloNuclei::Process(Int_t ntr){ Line 398  void CaloNuclei::Process(Int_t ntr){
398            postq += mip;            postq += mip;
399          } else {          } else {
400            preq += mip;            preq += mip;
401            if (  strip == (track->GetCaloTrack()->tibar[plane][view]-1) || strip == (track->GetCaloTrack()->tibar[plane][view]-2) || strip == (track->GetCaloTrack()->tibar[plane][view]) ){            if ( ntr >= 0 ){
402              if ( qsplane != plane || qsview != view ){              if (  strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) ){
403                qsplane = plane;                if ( qsplane != plane || qsview != view ){
404                qsview = view;                  qsplane = plane;
405                ind++;                  qsview = view;
406                if ( debug && ind > 199 ) printf(" AAAGH!! \n");                  ind++;
407                qme[ind] = 0.;                  if ( debug && ind > 199 ) printf(" AAAGH!! \n");
408                    qme[ind] = 0.;
409                  };
410                  qme[ind] += mip;
411              };              };
412              qme[ind] += mip;              for ( Int_t ns = 0; ns < R ; ns++){
413            };                Int_t ms =  track->tibar[plane][view] - 1 - ns + (R - 1)/2;
414            for ( Int_t ns = 0; ns < R ; ns++){                if ( strip == ms ){
415              Int_t ms =  track->GetCaloTrack()->tibar[plane][view] - 1 - ns + (R - 1)/2;                  if ( qsplane2 != plane || qsview2 != view ){
416              if ( strip == ms ){                    qsplane2 = plane;
417                if ( qsplane2 != plane || qsview2 != view ){                    qsview2 = view;
418                  qsplane2 = plane;                    ind2++;
419                  qsview2 = view;                    if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
420                  ind2++;                    qme2[ind2] = 0.;
421                  if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");                  };
422                  qme2[ind2] = 0.;                  qme2[ind2] += mip;
423                };                };
424                qme2[ind2] += mip;              };  
425              } else {
426                if (  strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) ){
427                  if ( qsplane != plane || qsview != view ){
428                    qsplane = plane;
429                    qsview = view;
430                    ind++;
431                    if ( debug && ind > 199 ) printf(" AAAGH!! \n");
432                    qme[ind] = 0.;
433                  };
434                  qme[ind] += mip;
435              };              };
436            };                  for ( Int_t ns = 0; ns < R ; ns++){
437                  Int_t ms =  L2->GetCaloLevel2()->cibar[plane][view] - 1 - ns + (R - 1)/2;
438                  if ( strip == ms ){
439                    if ( qsplane2 != plane || qsview2 != view ){
440                      qsplane2 = plane;
441                      qsview2 = view;
442                      ind2++;
443                      if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
444                      qme2[ind2] = 0.;
445                    };
446                    qme2[ind2] += mip;
447                  };
448                };  
449              };
450          };                };      
451          //          //
452        };        };
# Line 315  void CaloNuclei::Process(Int_t ntr){ Line 466  void CaloNuclei::Process(Int_t ntr){
466      Long64_t work[200];      Long64_t work[200];
467      ind = 0;      ind = 0;
468      Int_t l = 0;      Int_t l = 0;
469        Int_t RN = 0;
470      Float_t qm = 0.;      Float_t qm = 0.;
471      Float_t qm2 = 0.;      Float_t qm2 = 0.;
472      //      //
# Line 325  void CaloNuclei::Process(Int_t ntr){ Line 477  void CaloNuclei::Process(Int_t ntr){
477      while ( l < uplim && ind < interplane ){      while ( l < uplim && ind < interplane ){
478        qm = TMath::KOrdStat(interplane,qme,ind,work);        qm = TMath::KOrdStat(interplane,qme,ind,work);
479        if ( qm >= qmt ){        if ( qm >= qmt ){
480          if ( l < 3 ) qpremean += qm;          if ( l < 3 ){
481              qpremean += qm;
482              RN++;
483            };
484          l++;          l++;
485          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);
486        };        };
487        ind++;        ind++;
488      };      };
489      //      //
490      qpremean /= l;      qpremean /= (Float_t)RN;
491      //      //
492      ind = 0;      ind = 0;
493      l = 0;      l = 0;
494        RN = 0;
495      while ( l < uplim && ind < interplane ){      while ( l < uplim && ind < interplane ){
496        qm2 = TMath::KOrdStat(interplane,qme2,ind,work);        qm2 = TMath::KOrdStat(interplane,qme2,ind,work);
497        if ( qm2 >= qmt ){                if ( qm2 >= qmt ){        
498          if ( l < N ) qpremeanN += qm2;          if ( l < N ){
499              qpremeanN += qm2;
500              RN++;
501            };
502          l++;          l++;
503          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);
504        };        };
505        ind++;        ind++;
506      };      };
507      //      //
508      qpremeanN /= l;      qpremeanN /= (Float_t)RN;
509        UN = RN;
510      //      //
511      if ( debug ) printf(" charge is %f \n",sqrt(qpremean));      if ( debug ) printf(" charge is %f \n",sqrt(qpremean));
512      //      //
# Line 358  void CaloNuclei::Process(Int_t ntr){ Line 518  void CaloNuclei::Process(Int_t ntr){
518          mesethr = mesethr2;          mesethr = mesethr2;
519        };        };
520        aldone = true;        aldone = true;
521        if ( mesethr > defethr ) goto retry;        if ( mesethr > defethr ){
522            interplane = 0;
523            preq = 0.;
524            postq = 0.;
525            qpremean = 0.;
526            qpremeanN = 0.;
527            multhit = false;
528            gap = false;
529            goto retry;
530          };
531      };      };
532    };    };
533    //    //
534      if ( debug ) this->Print();
535    if ( debug ) printf(" esci \n");    if ( debug ) printf(" esci \n");
536    //    //
537  };  };
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.23