/[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.4 by mocchiut, Thu Apr 5 13:38:32 2007 UTC
# Line 20  CaloNuclei::CaloNuclei(PamLevel2 *l2p){ Line 20  CaloNuclei::CaloNuclei(PamLevel2 *l2p){
20    PKT = 0;    PKT = 0;
21    atime = 0;    atime = 0;
22    N = 5;    N = 5;
23      R = 3;
24      //
25      debug = false;
26    //    //
27  };  };
28    
29  void CaloNuclei::Clear(){  void CaloNuclei::Clear(){
30    //    //
31      tr = 0;
32    interplane = 0;    interplane = 0;
33    preq = 0.;    preq = 0.;
34    postq = 0.;    postq = 0.;
# Line 40  void CaloNuclei::Clear(){ Line 44  void CaloNuclei::Clear(){
44    
45  void CaloNuclei::Print(){  void CaloNuclei::Print(){
46    //    //
47    printf("==============================\n");    Process();
48    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);    //
49    printf(" interplane: %i\n",interplane);    printf("========================================================================\n");
50    printf(" ethr: %f \n",ethr);    printf(" OBT: %u PKT: %u ATIME: %u Track %i \n",OBT,PKT,atime,tr);
51    printf(" dedx1: %f \n",dedx1);    printf(" interplane [number of available dE/dx before interaction]: %i\n",interplane);
52    printf(" dedx3: %f \n",dedx3);    printf(" ethr [threshold used to determine interplane]:............ %f \n",ethr);
53    printf(" multhit: %i \n",multhit);    printf(" dedx1 [dE/dx from the first calorimeter plane]:........... %f \n",dedx1);
54    printf(" gap: %i \n",gap);    printf(" dedx3 [dE/dx (average) if the first 3 Si planes]:......... %f \n",dedx3);
55    printf(" preq: %f \n",preq);    printf(" multhit [true if interplane determined by multiple hits]:. %i \n",multhit);
56    printf(" postq: %f \n",postq);    printf(" gap [true if interplane determined by a gap]:............. %i \n",gap);
57    printf(" qpremean: %f \n",qpremean);    printf(" preq [total energy in MIP before the interaction plane]:.. %f \n",preq);
58    printf(" N: %i \n",N);    printf(" postq [total energy in MIP after the interaction plane]:.. %f \n",postq);
59    printf(" qpremeanN: %f \n",qpremeanN);    printf(" qpremean [truncated mean using 3 planes and 3 strips]:.... %f \n",qpremean);
60    printf("==============================\n");    printf(" N [no of used plane]:..................................... %i \n",N);
61      printf(" R [no strip used per plane ]:............................. %i \n",R);
62      printf(" qpremeanN [truncated mean using N planes and R strips]:... %f \n",qpremeanN);
63      printf("========================================================================\n");
64    //    //
65  };  };
66    
# Line 64  void CaloNuclei::Delete(){ Line 71  void CaloNuclei::Delete(){
71    
72    
73  void CaloNuclei::Process(){  void CaloNuclei::Process(){
74    //    Process(0);
75    };
76    
77    void CaloNuclei::Process(Int_t ntr){
78      //  
79    if ( !L2 ){    if ( !L2 ){
80      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");
81      printf(" ERROR: CaloNuclei variables not filled \n");      printf(" ERROR: CaloNuclei variables not filled \n");
# Line 86  void CaloNuclei::Process(){ Line 97  void CaloNuclei::Process(){
97    //    //
98    if ( !newentry ) return;    if ( !newentry ) return;
99    //    //
100    //  printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);    tr = ntr;
101      //
102      if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
103    //    //
104    Clear();    Clear();
105    //    //
106    PamTrack *track = 0;    PamTrack *track = 0;
107    track = L2->GetTrack(0);    track = L2->GetTrack(ntr);
108    //    //
109    Int_t view = 0;    Int_t view = 0;
110    Int_t plane = 0;    Int_t plane = 0;
# Line 133  void CaloNuclei::Process(){ Line 146  void CaloNuclei::Process(){
146    //    //
147   retry:   retry:
148    //    //
149    //  printf("retry\n");    if ( debug ) printf("retry\n");
150    //    //
151    interplane = 0;    interplane = 0;
152    //    //
# Line 157  void CaloNuclei::Process(){ Line 170  void CaloNuclei::Process(){
170      if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&      if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
171           ( strip == (track->GetCaloTrack()->tibar[plane][view]-1) || strip == (track->GetCaloTrack()->tibar[plane][view]-2) || strip == (track->GetCaloTrack()->tibar[plane][view]) )           ( strip == (track->GetCaloTrack()->tibar[plane][view]-1) || strip == (track->GetCaloTrack()->tibar[plane][view]-2) || strip == (track->GetCaloTrack()->tibar[plane][view]) )
172           && true ){                 && true ){      
173        //      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 ( 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]);
174        interpl[view] = plane;        interpl[view] = plane;
175        interv[view] = view;        interv[view] = view;
176        if ( splane[view] != plane || sview[view] != view ){        if ( splane[view] != plane || sview[view] != view ){
177          if ( nhit[view] > 1 ){          if ( nhit[view] > 1 ){
178            wmulthit[view] = true;            wmulthit[view] = true;
179              //      if ( splane[view] == -1 ) splane[view] = 0; //
180              //      if ( sview[view] == -1 ) sview[view] = view; //
181            interpl[view] = splane[view];            interpl[view] = splane[view];
182            interv[view] = sview[view];            interv[view] = sview[view];    
           //      ii = L2->GetCaloLevel1()->istrip;  
183          };          };
184          if ( plane > splane[view]+gapth ){          if ( plane > splane[view]+gapth ){
185            wgap[view] = true;            wgap[view] = true;
186              //      if ( splane[view] == -1 ) splane[view] = 0;//
187              //      if ( sview[view] == -1 ) sview[view] = view; //
188            interpl[view] = splane[view];            interpl[view] = splane[view];
189            interv[view] = sview[view];            interv[view] = sview[view];
           //      ii = L2->GetCaloLevel1()->istrip;  
190          };          };
191          splane[view] = plane;          splane[view] = plane;
192          sview[view] = view;          sview[view] = view;
# Line 185  void CaloNuclei::Process(){ Line 200  void CaloNuclei::Process(){
200      //          //    
201    };    };
202    //    //
203    //  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);
204    Int_t winterplane[2] = {-1,-1};    Int_t winterplane[2] = {-1,-1};
205    //    //
206    for ( Int_t view = 0; view < 2; view++){    for ( Int_t view = 0; view < 2; view++){
# Line 224  void CaloNuclei::Process(){ Line 239  void CaloNuclei::Process(){
239      };      };
240    };    };
241    //      //  
242    //  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);
243    //  printf("3conversion winterpl0 %i winterpl1 %i \n",winterplane[0],winterplane[1]);    if ( debug ) printf("3conversion winterpl0 %i winterpl1 %i \n",winterplane[0],winterplane[1]);
244    //    //
245    Int_t ipl = 0;    Int_t ipl = 0;
246    if ( interplane > 0 ){    if ( interplane > 0 ){
# Line 236  void CaloNuclei::Process(){ Line 251  void CaloNuclei::Process(){
251      Int_t ind = -1;      Int_t ind = -1;
252      Int_t qsplane = -1;      Int_t qsplane = -1;
253      Int_t qsview = -1;      Int_t qsview = -1;
254        Int_t ind2 = -1;
255        Int_t qsplane2 = -1;
256        Int_t qsview2 = -1;
257      Float_t qme[200];      Float_t qme[200];
258      memset(qme,0,200*sizeof(Float_t));      memset(qme,0,200*sizeof(Float_t));
259        Float_t qme2[2112];
260        memset(qme2,0,2112*sizeof(Float_t));
261      //      //
262      while ( ii<L2->GetCaloLevel1()->istrip ){      while ( ii<L2->GetCaloLevel1()->istrip ){
263        //        //
# Line 258  void CaloNuclei::Process(){ Line 278  void CaloNuclei::Process(){
278                qsplane = plane;                qsplane = plane;
279                qsview = view;                qsview = view;
280                ind++;                ind++;
281                if ( ind > 199 ) printf(" AAAGH!! \n");                if ( debug && ind > 199 ) printf(" AAAGH!! \n");
282                qme[ind] = 0.;                qme[ind] = 0.;
283              };              };
284              qme[ind] += mip;              qme[ind] += mip;
285            };            };
286              for ( Int_t ns = 0; ns < R ; ns++){
287                Int_t ms =  track->GetCaloTrack()->tibar[plane][view] - 1 - ns + (R - 1)/2;
288                if ( strip == ms ){
289                  if ( qsplane2 != plane || qsview2 != view ){
290                    qsplane2 = plane;
291                    qsview2 = view;
292                    ind2++;
293                    if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
294                    qme2[ind2] = 0.;
295                  };
296                  qme2[ind2] += mip;
297                };
298              };    
299          };                };      
300          //          //
301        };        };
# Line 273  void CaloNuclei::Process(){ Line 306  void CaloNuclei::Process(){
306      //      //
307      // 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...
308      //          //    
309      //    for (Int_t l=0; l < interplane; l++){      if ( debug ){
310      //       printf(" qme[%i] = %f \n",l,qme[l]);          for (Int_t l=0; l < interplane; l++){
311      //       //      if ( qme[ind] < 6.25 ) qme[ind] = 10000.;                printf(" qme[%i] = %f \n",l,qme[l]);  
312      //     };        };
313        };
314        //
315      Long64_t work[200];      Long64_t work[200];
316      ind = 0;      ind = 0;
317      Int_t l = 0;      Int_t l = 0;
318      Float_t qm = 0.;      Float_t qm = 0.;
319        Float_t qm2 = 0.;
320        //
321        Float_t qmt = ethr*0.8; // *0.9
322        //
323      Int_t uplim = TMath::Max(3,N);      Int_t uplim = TMath::Max(3,N);
324        //
325      while ( l < uplim && ind < interplane ){      while ( l < uplim && ind < interplane ){
326        qm = TMath::KOrdStat(interplane,qme,ind,work);        qm = TMath::KOrdStat(interplane,qme,ind,work);
327        if ( qm >= mesethr*0.9 ){        if ( qm >= qmt ){
328          if ( l < 3 ) qpremean += qm;          if ( l < 3 ) qpremean += qm;
         if ( l < N ) qpremeanN += qm;  
329          l++;          l++;
330          //      printf(" value no %i qm %f \n",l,qm);          if ( debug ) printf(" value no %i qm %f qmt %f \n",l,qm,qmt);
331        };        };
332        ind++;        ind++;
333      };      };
334      //      //
335      qpremean /= l;      qpremean /= l;
     qpremeanN /= N;  
336      //      //
337      //  printf(" charge is %f \n",sqrt(qpremean));      ind = 0;
338        l = 0;
339        while ( l < uplim && ind < interplane ){
340          qm2 = TMath::KOrdStat(interplane,qme2,ind,work);
341          if ( qm2 >= qmt ){        
342            if ( l < N ) qpremeanN += qm2;
343            l++;
344            if ( debug ) printf(" qm2 value no %i qm %f qmt %f \n",l,qm2,qmt);
345          };
346          ind++;
347        };
348        //
349        qpremeanN /= l;
350      //      //
351      //    printf("preq postq\n");          if ( debug ) printf(" charge is %f \n",sqrt(qpremean));
352      //      //
353      if ( mesethr != ethr && interplane >= 3 && !aldone ){      if ( mesethr != ethr && interplane >= 3 && !aldone ){
354        Float_t mesethr2 = (sqrt(qpremean) - 0.50)*(sqrt(qpremean) - 0.50);        Float_t mesethr2 = (sqrt(qpremean) - 0.50)*(sqrt(qpremean) - 0.50);
# Line 312  void CaloNuclei::Process(){ Line 362  void CaloNuclei::Process(){
362      };      };
363    };    };
364    //    //
365    //  printf(" esci \n");    if ( debug ) printf(" esci \n");
366    //    //
367  };  };
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.23