/[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.7 by mocchiut, Thu May 31 10:10:52 2007 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    tr = 0;    tr = 0;
33      sntr = 0;
34    interplane = 0;    interplane = 0;
35    preq = 0.;    preq = 0.;
36    postq = 0.;    postq = 0.;
# Line 47  void CaloNuclei::Print(){ Line 49  void CaloNuclei::Print(){
49    Process();    Process();
50    //    //
51    printf("========================================================================\n");    printf("========================================================================\n");
52    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);
53    printf(" interplane [number of available dE/dx before interaction]: %i\n",interplane);    printf(" interplane [number of available dE/dx before interaction]:.. %i\n",interplane);
54    printf(" ethr [threshold used to determine interplane]:............ %f \n",ethr);    printf(" ethr [threshold used to determine interplane]:.............. %f \n",ethr);
55    printf(" dedx1 [dE/dx from the first calorimeter plane]:........... %f \n",dedx1);    printf(" dedx1 [dE/dx from the first calorimeter plane]:............. %f \n",dedx1);
56    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);
57    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);
58    printf(" gap [true if interplane determined by a gap]:............. %i \n",gap);    printf(" multhit [true if interplane determined by multiple hits]:... %i \n",multhit);
59    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);
60    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);
61    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);
62    printf(" N [no of used plane]:..................................... %i \n",N);    printf(" qpremean [truncated mean using 3 planes and 3 strips]:...... %f \n",qpremean);
63    printf(" R [no strip used per plane ]:............................. %i \n",R);    printf(" N [no of used plane]:....................................... %i \n",N);
64    printf(" qpremeanN [truncated mean using N planes and R strips]:... %f \n",qpremeanN);    printf(" R [no strip used per plane ]:............................... %i \n",R);
65      printf(" qpremeanN [truncated mean using N planes and R strips]:..... %f \n",qpremeanN);
66    printf("========================================================================\n");    printf("========================================================================\n");
67    //    //
68  };  };
# Line 85  void CaloNuclei::Process(Int_t ntr){ Line 88  void CaloNuclei::Process(Int_t ntr){
88    Bool_t newentry = false;    Bool_t newentry = false;
89    //    //
90    if ( L2->IsORB() ){    if ( L2->IsORB() ){
91      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 ){
92        newentry = true;        newentry = true;
93        OBT = L2->GetOrbitalInfo()->OBT;        OBT = L2->GetOrbitalInfo()->OBT;
94        PKT = L2->GetOrbitalInfo()->pkt_num;        PKT = L2->GetOrbitalInfo()->pkt_num;
95        atime = L2->GetOrbitalInfo()->absTime;        atime = L2->GetOrbitalInfo()->absTime;
96          sntr = ntr;
97      };      };
98    } else {    } else {
99      newentry = true;      newentry = true;
# Line 102  void CaloNuclei::Process(Int_t ntr){ Line 106  void CaloNuclei::Process(Int_t ntr){
106    if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);    if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
107    //    //
108    Clear();    Clear();
109      if ( debug ) printf(" Always calculate stdedx1 \n");
110    //    //
111    PamTrack *track = 0;    // Always calculate stdedx1
   track = L2->GetTrack(ntr);  
112    //    //
113    Int_t view = 0;    Int_t view = 0;
114    Int_t plane = 0;    Int_t plane = 0;
115    Int_t strip = 0;    Int_t strip = 0;
116      Int_t indx = 0;
117      Float_t vfpl[96];
118      Int_t stfpl[96];
119      memset(vfpl, 0, 96*sizeof(Float_t));
120      memset(stfpl, 0, 96*sizeof(Int_t));
121    Float_t mip = 0.;    Float_t mip = 0.;
122      for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
123        //
124        mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
125        //
126        // put in vfpl vector the energy release on the first plane
127        //
128        if ( strip != -1 && view == 1 && plane == 0 ) {
129          stfpl[indx] = strip;
130          vfpl[indx] = mip;
131          indx++;
132        };
133        //
134      };
135      //
136      if ( debug ) printf(" find energy released along the strip of maximum on the first plane and on the two neighbour strips \n");
137      //
138      // find energy released along the strip of maximum on the first plane and on the two neighbour strips
139      //
140      if ( indx > 0 ){
141        Int_t mindx = (Int_t)TMath::LocMax(indx,stfpl);
142        for (Int_t ii=0; ii<indx; ii++){
143          if ( stfpl[ii] == stfpl[mindx] ) stdedx1 += vfpl[ii];
144          if ( (mindx-1)>=0 && stfpl[ii] == stfpl[mindx-1] ) stdedx1 += vfpl[ii];
145          if ( (mindx+1)<96 && stfpl[ii] == stfpl[mindx+1] ) stdedx1 += vfpl[ii];
146        };
147      } else {
148        stdedx1 = 0.;
149      };
150      //
151      if ( debug ) printf(" if ( !usetrack ) return: usetrack %i ntr %i \n",usetrack,ntr);
152      //
153      //
154      if ( !usetrack ) return;
155      //
156      PamTrack *ptrack = 0;
157      CaloTrkVar *track = 0;
158      if ( ntr >= 0 ){
159        ptrack = L2->GetTrack(ntr);
160        if ( ptrack ) track = ptrack->GetCaloTrack();
161      } else {
162        track = L2->GetCaloStoredTrack(ntr);
163      };
164      //
165      if ( !track && ntr >= 0 ){
166        printf(" ERROR: cannot find any track!\n");
167        printf(" ERROR: CaloNuclei variables not completely filled \n");
168        return;  
169      };
170      //
171    //  Float_t defethr = 6. * 0.90;    //  Float_t defethr = 6. * 0.90;
172    Float_t defethr = 6.25; // = (sqrt(9) - 0.5) ** 2.;    Float_t defethr = 6.25; // = (sqrt(9) - 0.5) ** 2.;
173    //    //
# Line 119  void CaloNuclei::Process(Int_t ntr){ Line 177  void CaloNuclei::Process(Int_t ntr){
177      //      //
178      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
179      //      //
180      if ( strip != -1 &&      if ( ntr >= 0 ){
181           view == 1   &&        //
182           plane == 0  &&        if ( strip != -1 &&
183           ( strip == (track->GetCaloTrack()->tibar[0][1]-1) || strip == (track->GetCaloTrack()->tibar[0][1]-2) || strip == (track->GetCaloTrack()->tibar[0][1]) )             view == 1   &&
184           && true ){                   plane == 0  &&
185        dedx1 += mip;             ( strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) || strip == (track->tibar[0][1]) )
186      };             && true ){      
187      if ( strip != -1 &&          dedx1 += mip;
188           (( view == 1 && ( plane == 0 || plane == 1 ) ) ||        };
189           ( view == 0 && plane == 0 )) &&        if ( strip != -1 &&
190           (( 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 ) ) ||
191            ( 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 )) &&
192            ( 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) )) ||
193           true ){              ( view == 1 && ( strip == track->tibar[0][1] || strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) )) ||
194        dedx3 += mip;              ( view == 1 && ( strip == track->tibar[1][1] || strip == (track->tibar[1][1]-1) || strip == (track->tibar[1][1]-2) ))) &&
195               true ){
196            dedx3 += mip;
197          };
198        } else {
199          //
200          if ( strip != -1 &&
201               view == 1   &&
202               plane == 0  &&
203               ( strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) || strip == (L2->GetCaloLevel2()->cibar[0][1]) )
204               && true ){      
205            dedx1 += mip;
206          };
207          if ( strip != -1 &&
208               (( view == 1 && ( plane == 0 || plane == 1 ) ) ||
209                ( view == 0 && plane == 0 )) &&
210               (( view == 0 && ( strip == L2->GetCaloLevel2()->cibar[0][0] || strip == (L2->GetCaloLevel2()->cibar[0][0]-1) || strip == (L2->GetCaloLevel2()->cibar[0][0]-2) )) ||
211                ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[0][1] || strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) )) ||
212                ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[1][1] || strip == (L2->GetCaloLevel2()->cibar[1][1]-1) || strip == (L2->GetCaloLevel2()->cibar[1][1]-2) ))) &&
213               true ){
214            dedx3 += mip;
215          };
216      };      };
217      //          //    
218    };    };
# Line 167  void CaloNuclei::Process(Int_t ntr){ Line 246  void CaloNuclei::Process(Int_t ntr){
246      //      //
247      mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);          mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);    
248      //      //
249      if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&      if ( ntr >= 0 ){
250           ( 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] &&
251           && true ){                   ( strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) )
252        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 ){      
253        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]);
254        interv[view] = view;          interpl[view] = plane;
255        if ( splane[view] != plane || sview[view] != view ){          interv[view] = view;
256          if ( nhit[view] > 1 ){          if ( splane[view] != plane || sview[view] != view ){
257            wmulthit[view] = true;            if ( nhit[view] > 1 ){
258            //      if ( splane[view] == -1 ) splane[view] = 0; //              wmulthit[view] = true;
259            //      if ( sview[view] == -1 ) sview[view] = view; //              //    if ( splane[view] == -1 ) splane[view] = 0; //
260            interpl[view] = splane[view];              //    if ( sview[view] == -1 ) sview[view] = view; //
261            interv[view] = sview[view];                  interpl[view] = splane[view];
262                interv[view] = sview[view];  
263          };          };
264          if ( plane > splane[view]+gapth ){            if ( plane > splane[view]+gapth ){
265            wgap[view] = true;              wgap[view] = true;
266            //      if ( splane[view] == -1 ) splane[view] = 0;//              //    if ( splane[view] == -1 ) splane[view] = 0;//
267            //      if ( sview[view] == -1 ) sview[view] = view; //              //    if ( sview[view] == -1 ) sview[view] = view; //
268            interpl[view] = splane[view];              interpl[view] = splane[view];
269            interv[view] = sview[view];              interv[view] = sview[view];
270              };
271              splane[view] = plane;
272              sview[view] = view;
273              nhit[view] = 1;
274            } else {
275              nhit[view]++;
276            };
277          };
278        } else {
279          if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
280               ( strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) )
281               && true ){      
282            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]);
283            interpl[view] = plane;
284            interv[view] = view;
285            if ( splane[view] != plane || sview[view] != view ){
286              if ( nhit[view] > 1 ){
287                wmulthit[view] = true;
288                //    if ( splane[view] == -1 ) splane[view] = 0; //
289                //    if ( sview[view] == -1 ) sview[view] = view; //
290                interpl[view] = splane[view];
291                interv[view] = sview[view];  
292            };
293              if ( plane > splane[view]+gapth ){
294                wgap[view] = true;
295                //    if ( splane[view] == -1 ) splane[view] = 0;//
296                //    if ( sview[view] == -1 ) sview[view] = view; //
297                interpl[view] = splane[view];
298                interv[view] = sview[view];
299              };
300              splane[view] = plane;
301              sview[view] = view;
302              nhit[view] = 1;
303            } else {
304              nhit[view]++;
305          };          };
         splane[view] = plane;  
         sview[view] = view;  
         nhit[view] = 1;  
       } else {  
         nhit[view]++;  
306        };        };
307      };      };
308      //      //
# Line 273  void CaloNuclei::Process(Int_t ntr){ Line 383  void CaloNuclei::Process(Int_t ntr){
383            postq += mip;            postq += mip;
384          } else {          } else {
385            preq += mip;            preq += mip;
386            if (  strip == (track->GetCaloTrack()->tibar[plane][view]-1) || strip == (track->GetCaloTrack()->tibar[plane][view]-2) || strip == (track->GetCaloTrack()->tibar[plane][view]) ){            if ( ntr >= 0 ){
387              if ( qsplane != plane || qsview != view ){              if (  strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) ){
388                qsplane = plane;                if ( qsplane != plane || qsview != view ){
389                qsview = view;                  qsplane = plane;
390                ind++;                  qsview = view;
391                if ( debug && ind > 199 ) printf(" AAAGH!! \n");                  ind++;
392                qme[ind] = 0.;                  if ( debug && ind > 199 ) printf(" AAAGH!! \n");
393                    qme[ind] = 0.;
394                  };
395                  qme[ind] += mip;
396              };              };
397              qme[ind] += mip;              for ( Int_t ns = 0; ns < R ; ns++){
398            };                Int_t ms =  track->tibar[plane][view] - 1 - ns + (R - 1)/2;
399            for ( Int_t ns = 0; ns < R ; ns++){                if ( strip == ms ){
400              Int_t ms =  track->GetCaloTrack()->tibar[plane][view] - 1 - ns + (R - 1)/2;                  if ( qsplane2 != plane || qsview2 != view ){
401              if ( strip == ms ){                    qsplane2 = plane;
402                if ( qsplane2 != plane || qsview2 != view ){                    qsview2 = view;
403                  qsplane2 = plane;                    ind2++;
404                  qsview2 = view;                    if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
405                  ind2++;                    qme2[ind2] = 0.;
406                  if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");                  };
407                  qme2[ind2] = 0.;                  qme2[ind2] += mip;
408                  };
409                };  
410              } else {
411                if (  strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) ){
412                  if ( qsplane != plane || qsview != view ){
413                    qsplane = plane;
414                    qsview = view;
415                    ind++;
416                    if ( debug && ind > 199 ) printf(" AAAGH!! \n");
417                    qme[ind] = 0.;
418                };                };
419                qme2[ind2] += mip;                qme[ind] += mip;
420              };              };
421            };                  for ( Int_t ns = 0; ns < R ; ns++){
422                  Int_t ms =  L2->GetCaloLevel2()->cibar[plane][view] - 1 - ns + (R - 1)/2;
423                  if ( strip == ms ){
424                    if ( qsplane2 != plane || qsview2 != view ){
425                      qsplane2 = plane;
426                      qsview2 = view;
427                      ind2++;
428                      if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
429                      qme2[ind2] = 0.;
430                    };
431                    qme2[ind2] += mip;
432                  };
433                };  
434              };
435          };                };      
436          //          //
437        };        };
# Line 315  void CaloNuclei::Process(Int_t ntr){ Line 451  void CaloNuclei::Process(Int_t ntr){
451      Long64_t work[200];      Long64_t work[200];
452      ind = 0;      ind = 0;
453      Int_t l = 0;      Int_t l = 0;
454        Int_t RN = 0;
455      Float_t qm = 0.;      Float_t qm = 0.;
456      Float_t qm2 = 0.;      Float_t qm2 = 0.;
457      //      //
# Line 325  void CaloNuclei::Process(Int_t ntr){ Line 462  void CaloNuclei::Process(Int_t ntr){
462      while ( l < uplim && ind < interplane ){      while ( l < uplim && ind < interplane ){
463        qm = TMath::KOrdStat(interplane,qme,ind,work);        qm = TMath::KOrdStat(interplane,qme,ind,work);
464        if ( qm >= qmt ){        if ( qm >= qmt ){
465          if ( l < 3 ) qpremean += qm;          if ( l < 3 ){
466              qpremean += qm;
467              RN++;
468            };
469          l++;          l++;
470          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);
471        };        };
472        ind++;        ind++;
473      };      };
474      //      //
475      qpremean /= l;      qpremean /= (Float_t)RN;
476      //      //
477      ind = 0;      ind = 0;
478      l = 0;      l = 0;
479        RN = 0;
480      while ( l < uplim && ind < interplane ){      while ( l < uplim && ind < interplane ){
481        qm2 = TMath::KOrdStat(interplane,qme2,ind,work);        qm2 = TMath::KOrdStat(interplane,qme2,ind,work);
482        if ( qm2 >= qmt ){                if ( qm2 >= qmt ){        
483          if ( l < N ) qpremeanN += qm2;          if ( l < N ){
484              qpremeanN += qm2;
485              RN++;
486            };
487          l++;          l++;
488          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);
489        };        };
490        ind++;        ind++;
491      };      };
492      //      //
493      qpremeanN /= l;      qpremeanN /= (Float_t)RN;
494      //      //
495      if ( debug ) printf(" charge is %f \n",sqrt(qpremean));      if ( debug ) printf(" charge is %f \n",sqrt(qpremean));
496      //      //
# Line 358  void CaloNuclei::Process(Int_t ntr){ Line 502  void CaloNuclei::Process(Int_t ntr){
502          mesethr = mesethr2;          mesethr = mesethr2;
503        };        };
504        aldone = true;        aldone = true;
505        if ( mesethr > defethr ) goto retry;        if ( mesethr > defethr ){
506            interplane = 0;
507            preq = 0.;
508            postq = 0.;
509            qpremean = 0.;
510            qpremeanN = 0.;
511            multhit = false;
512            gap = false;
513            goto retry;
514          };
515      };      };
516    };    };
517    //    //
518      if ( debug ) this->Print();
519    if ( debug ) printf(" esci \n");    if ( debug ) printf(" esci \n");
520    //    //
521  };  };
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.23