CaloNuclei/src/CaloNuclei.cpp

#include <CaloNuclei.h>

//--------------------------------------
/**
 * Default constructor 
 */
CaloNuclei::CaloNuclei(){
  Clear();
};

CaloNuclei::CaloNuclei(PamLevel2 *l2p){  
  //
  Clear();
  //
  L2 = l2p;
  //
  if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
  //
  OBT = 0;
  PKT = 0;
  atime = 0;
  N = 5;
  //
};

void CaloNuclei::Clear(){
  //
  interplane = 0;
  preq = 0.;
  postq = 0.;
  dedx1 = 0.;
  dedx3 = 0.;
  qpremean = 0.; 
  qpremeanN = 0.; 
  //
  multhit = false;
  gap = false;
  //
};

void CaloNuclei::Print(){
  //
  printf("==============================\n");
  printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
  printf(" interplane: %i\n",interplane);
  printf(" ethr: %f \n",ethr); 
  printf(" dedx1: %f \n",dedx1);
  printf(" dedx3: %f \n",dedx3);
  printf(" multhit: %i \n",multhit);
  printf(" gap: %i \n",gap);
  printf(" preq: %f \n",preq);
  printf(" postq: %f \n",postq);
  printf(" qpremean: %f \n",qpremean); 
  printf(" N: %i \n",N); 
  printf(" qpremeanN: %f \n",qpremeanN); 
  printf("==============================\n");
  //
};

void CaloNuclei::Delete(){
  Clear();
  //delete this;
};


void CaloNuclei::Process(){
  //
  if ( !L2 ){
    printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
    printf(" ERROR: CaloNuclei variables not filled \n");
    return;
  };
  //
  Bool_t newentry = false;
  //
  if ( L2->IsORB() ){
    if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
      newentry = true;
      OBT = L2->GetOrbitalInfo()->OBT;
      PKT = L2->GetOrbitalInfo()->pkt_num;
      atime = L2->GetOrbitalInfo()->absTime;
    };
  } else {
    newentry = true;
  };
  //
  if ( !newentry ) return;
  //
  //  printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
  //
  Clear();
  //
  PamTrack *track = 0;
  track = L2->GetTrack(0);
  //
  Int_t view = 0;
  Int_t plane = 0;
  Int_t strip = 0;
  Float_t mip = 0.;
  //  Float_t defethr = 6. * 0.90;
  Float_t defethr = 6.25; // = (sqrt(9) - 0.5) ** 2.;
  //
  // Calculate dedx1 and dedx3
  //
  for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
    //
    mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
    //
    if ( strip != -1 && 
         view == 1   && 
         plane == 0  &&
         ( strip == (track->GetCaloTrack()->tibar[0][1]-1) || strip == (track->GetCaloTrack()->tibar[0][1]-2) || strip == (track->GetCaloTrack()->tibar[0][1]) ) 
         && true ){      
      dedx1 += mip;
    };
    if ( strip != -1 && 
         (( view == 1 && ( plane == 0 || plane == 1 ) ) ||
         ( view == 0 && plane == 0 )) &&
         (( view == 0 && ( strip == track->GetCaloTrack()->tibar[0][0] || strip == (track->GetCaloTrack()->tibar[0][0]-1) || strip == (track->GetCaloTrack()->tibar[0][0]-2) )) ||
          ( view == 1 && ( strip == track->GetCaloTrack()->tibar[0][1] || strip == (track->GetCaloTrack()->tibar[0][1]-1) || strip == (track->GetCaloTrack()->tibar[0][1]-2) )) ||
          ( view == 1 && ( strip == track->GetCaloTrack()->tibar[1][1] || strip == (track->GetCaloTrack()->tibar[1][1]-1) || strip == (track->GetCaloTrack()->tibar[1][1]-2) ))) &&
         true ){
      dedx3 += mip;
    };
    //    
  };
  //
  dedx3 /= 3.;
  //  Float_t mesethr = dedx1 * 0.90;
  Float_t mesethr = 0.;
  if ( dedx1 > 0. ) mesethr = (sqrt(dedx1) - 0.50)*(sqrt(dedx1) - 0.50);
  Bool_t aldone = false;
  //
 retry:
  // 
  //  printf("retry\n");
  //
  interplane = 0;
  //
  ethr = TMath::Max(defethr,mesethr);
  //
  // Find the interaction plane "interplane"
  //
  Int_t gapth = 3;
  Int_t nhit[2] = {0,0};
  Int_t splane[2] = {-1,-1};
  Int_t sview[2] = {-1,-1};
  Int_t interpl[2] = {-1,-1};
  Int_t interv[2] = {-1,-1};
  Bool_t wmulthit[2] = {false,false};
  Bool_t wgap[2] = {false,false};
  Int_t ii = 0;
  while ( ii<L2->GetCaloLevel1()->istrip ){
    //
    mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);    
    //
    if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
         ( strip == (track->GetCaloTrack()->tibar[plane][view]-1) || strip == (track->GetCaloTrack()->tibar[plane][view]-2) || strip == (track->GetCaloTrack()->tibar[plane][view]) ) 
         && true ){      
      //      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]);
      interpl[view] = plane;
      interv[view] = view;
      if ( splane[view] != plane || sview[view] != view ){
        if ( nhit[view] > 1 ){
          wmulthit[view] = true;
          interpl[view] = splane[view];
          interv[view] = sview[view];
          //      ii = L2->GetCaloLevel1()->istrip;
        };
        if ( plane > splane[view]+gapth ){
          wgap[view] = true;
          interpl[view] = splane[view];
          interv[view] = sview[view];
          //      ii = L2->GetCaloLevel1()->istrip;
        };
        splane[view] = plane;
        sview[view] = view;
        nhit[view] = 1;
      } else {
        nhit[view]++;
      };
    };
    //
    ii++;
    //    
  };
  //
  //  printf("conversion interpl %i interv %i multhit %i interplane %i \n",interpl[0],interv[0],multhit,interplane);
  Int_t winterplane[2] = {-1,-1};
  //
  for ( Int_t view = 0; view < 2; view++){
    //
    if ( nhit[view] > 1 && !wmulthit[view] && !wgap[view] ){
      wmulthit[view] = true;
      interpl[view] = splane[view];
      interv[view] = sview[view];
    };
    //
    if ( wmulthit[view] ) multhit = true;
    if ( wgap[view] ) gap = true;
    //
    // convert view and plane number of interaction plane into number of available dE/dx measurements before the interaction plane
    //    
    if ( interpl[view] >= 0 ) {
      if ( interv[view] == 0 ){
        winterplane[view] = (1 + interpl[view]) * 2;
      } else {
        winterplane[view] = (1 + interpl[view]) + (1 + interpl[view] - 1);
      };
      if ( wmulthit[view] ) winterplane[view]--;
    };
  };
  if ( winterplane[0] > 0 && winterplane[1] > 0 ){
    if ( multhit ){
      interplane = TMath::Min(winterplane[0],winterplane[1]);
    } else {
      interplane = TMath::Max(winterplane[0],winterplane[1]);
    };
  } else {
    if ( !winterplane[0] || !winterplane[1] ){
      interplane = 0;
    } else {
      interplane = TMath::Max(winterplane[0],winterplane[1]);
    };
  };
  //  
  //  printf("2conversion interpl %i interv %i multhit %i interplane %i \n",interpl[1],interv[1],multhit,interplane);
  //  printf("3conversion winterpl0 %i winterpl1 %i \n",winterplane[0],winterplane[1]);
  //
  Int_t ipl = 0;
  if ( interplane > 0 ){
    //
    // Calculate preq, postq, qpremean
    //
    ii = 0;
    Int_t ind = -1;
    Int_t qsplane = -1;
    Int_t qsview = -1;
    Float_t qme[200];
    memset(qme,0,200*sizeof(Float_t));
    //
    while ( ii<L2->GetCaloLevel1()->istrip ){
      //
      mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);    
      //
      if ( strip != -1 ){
        if ( view == 0 ){
          ipl = (1 + plane) * 2;
        } else {
          ipl = (1 + plane) + (1 + plane - 1 );
        };
        if ( ipl > interplane ){
          postq += mip;
        } else {
          preq += mip;
          if (  strip == (track->GetCaloTrack()->tibar[plane][view]-1) || strip == (track->GetCaloTrack()->tibar[plane][view]-2) || strip == (track->GetCaloTrack()->tibar[plane][view]) ){
            if ( qsplane != plane || qsview != view ){
              qsplane = plane;
              qsview = view;
              ind++;
              if ( ind > 199 ) printf(" AAAGH!! \n");
              qme[ind] = 0.;
            };
            qme[ind] += mip;
          };
        };      
        //
      };
      //
      ii++;
      //    
    };
    //
    // here we must calculate qpremean, order vector qme, select 3 lowest measurements and caculate the mean...
    //    
    //    for (Int_t l=0; l < interplane; l++){
    //       printf(" qme[%i] = %f \n",l,qme[l]);  
    //       //      if ( qme[ind] < 6.25 ) qme[ind] = 10000.;      
    //     };
    Long64_t work[200];
    ind = 0;
    Int_t l = 0;
    Float_t qm = 0.;
    Int_t uplim = TMath::Max(3,N);
    while ( l < uplim && ind < interplane ){
      qm = TMath::KOrdStat(interplane,qme,ind,work);
      if ( qm >= mesethr*0.9 ){ 
        if ( l < 3 ) qpremean += qm;
        if ( l < N ) qpremeanN += qm;
        l++;
        //      printf(" value no %i qm %f \n",l,qm); 
      };
      ind++;
    };
    //
    qpremean /= l;
    qpremeanN /= N;
    //
    //  printf(" charge is %f \n",sqrt(qpremean));
    //
    //    printf("preq postq\n");    
    //
    if ( mesethr != ethr && interplane >= 3 && !aldone ){
      Float_t mesethr2 = (sqrt(qpremean) - 0.50)*(sqrt(qpremean) - 0.50);
      if ( mesethr2 < mesethr*0.90 ){
        mesethr = (sqrt(dedx1) - 0.25)*(sqrt(dedx1) - 0.25);
      } else {
        mesethr = mesethr2;
      };
      aldone = true;
      if ( mesethr > defethr ) goto retry;
    };
  };
  //
  //  printf(" esci \n");
  //
};
1	mocchiut	1.1	#include <CaloNuclei.h>
2
3			//--------------------------------------
4			/**
5			* Default constructor
6			*/
7			CaloNuclei::CaloNuclei(){
8			Clear();
9			};
10
11			CaloNuclei::CaloNuclei(PamLevel2 *l2p){
12			//
13			Clear();
14			//
15			L2 = l2p;
16			//
17			if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
18			//
19			OBT = 0;
20			PKT = 0;
21			atime = 0;
22			N = 5;
23			//
24			};
25
26			void CaloNuclei::Clear(){
27			//
28			interplane = 0;
29			preq = 0.;
30			postq = 0.;
31			dedx1 = 0.;
32			dedx3 = 0.;
33			qpremean = 0.;
34			qpremeanN = 0.;
35			//
36			multhit = false;
37			gap = false;
38			//
39			};
40
41			void CaloNuclei::Print(){
42			//
43			printf("==============================\n");
44			printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
45			printf(" interplane: %i\n",interplane);
46			printf(" ethr: %f \n",ethr);
47			printf(" dedx1: %f \n",dedx1);
48			printf(" dedx3: %f \n",dedx3);
49			printf(" multhit: %i \n",multhit);
50			printf(" gap: %i \n",gap);
51			printf(" preq: %f \n",preq);
52			printf(" postq: %f \n",postq);
53			printf(" qpremean: %f \n",qpremean);
54			printf(" N: %i \n",N);
55			printf(" qpremeanN: %f \n",qpremeanN);
56			printf("==============================\n");
57			//
58			};
59
60			void CaloNuclei::Delete(){
61			Clear();
62			//delete this;
63			};
64
65
66			void CaloNuclei::Process(){
67			//
68			if ( !L2 ){
69			printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
70			printf(" ERROR: CaloNuclei variables not filled \n");
71			return;
72			};
73			//
74			Bool_t newentry = false;
75			//
76			if ( L2->IsORB() ){
77			if ( L2->GetOrbitalInfo()->pkt_num != PKT \|\| L2->GetOrbitalInfo()->OBT != OBT \|\| L2->GetOrbitalInfo()->absTime != atime ){
78			newentry = true;
79			OBT = L2->GetOrbitalInfo()->OBT;
80			PKT = L2->GetOrbitalInfo()->pkt_num;
81			atime = L2->GetOrbitalInfo()->absTime;
82			};
83			} else {
84			newentry = true;
85			};
86			//
87			if ( !newentry ) return;
88			//
89			// printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
90			//
91			Clear();
92			//
93			PamTrack *track = 0;
94			track = L2->GetTrack(0);
95			//
96			Int_t view = 0;
97			Int_t plane = 0;
98			Int_t strip = 0;
99			Float_t mip = 0.;
100			// Float_t defethr = 6. * 0.90;
101			Float_t defethr = 6.25; // = (sqrt(9) - 0.5) ** 2.;
102			//
103			// Calculate dedx1 and dedx3
104			//
105			for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
106			//
107			mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
108			//
109			if ( strip != -1 &&
110			view == 1 &&
111			plane == 0 &&
112			( strip == (track->GetCaloTrack()->tibar[0][1]-1) \|\| strip == (track->GetCaloTrack()->tibar[0][1]-2) \|\| strip == (track->GetCaloTrack()->tibar[0][1]) )
113			&& true ){
114			dedx1 += mip;
115			};
116			if ( strip != -1 &&
117			(( view == 1 && ( plane == 0 \|\| plane == 1 ) ) \|\|
118			( view == 0 && plane == 0 )) &&
119			(( view == 0 && ( strip == track->GetCaloTrack()->tibar[0][0] \|\| strip == (track->GetCaloTrack()->tibar[0][0]-1) \|\| strip == (track->GetCaloTrack()->tibar[0][0]-2) )) \|\|
120			( view == 1 && ( strip == track->GetCaloTrack()->tibar[0][1] \|\| strip == (track->GetCaloTrack()->tibar[0][1]-1) \|\| strip == (track->GetCaloTrack()->tibar[0][1]-2) )) \|\|
121			( view == 1 && ( strip == track->GetCaloTrack()->tibar[1][1] \|\| strip == (track->GetCaloTrack()->tibar[1][1]-1) \|\| strip == (track->GetCaloTrack()->tibar[1][1]-2) ))) &&
122			true ){
123			dedx3 += mip;
124			};
125			//
126			};
127			//
128			dedx3 /= 3.;
129			// Float_t mesethr = dedx1 * 0.90;
130			Float_t mesethr = 0.;
131			if ( dedx1 > 0. ) mesethr = (sqrt(dedx1) - 0.50)*(sqrt(dedx1) - 0.50);
132			Bool_t aldone = false;
133			//
134			retry:
135			//
136			// printf("retry\n");
137			//
138			interplane = 0;
139			//
140			ethr = TMath::Max(defethr,mesethr);
141			//
142			// Find the interaction plane "interplane"
143			//
144			Int_t gapth = 3;
145			Int_t nhit[2] = {0,0};
146			Int_t splane[2] = {-1,-1};
147			Int_t sview[2] = {-1,-1};
148			Int_t interpl[2] = {-1,-1};
149			Int_t interv[2] = {-1,-1};
150			Bool_t wmulthit[2] = {false,false};
151			Bool_t wgap[2] = {false,false};
152			Int_t ii = 0;
153			while ( ii<L2->GetCaloLevel1()->istrip ){
154			//
155			mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);
156			//
157			if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
158			( strip == (track->GetCaloTrack()->tibar[plane][view]-1) \|\| strip == (track->GetCaloTrack()->tibar[plane][view]-2) \|\| strip == (track->GetCaloTrack()->tibar[plane][view]) )
159			&& true ){
160			// 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]);
161			interpl[view] = plane;
162			interv[view] = view;
163			if ( splane[view] != plane \|\| sview[view] != view ){
164			if ( nhit[view] > 1 ){
165			wmulthit[view] = true;
166			interpl[view] = splane[view];
167			interv[view] = sview[view];
168			// ii = L2->GetCaloLevel1()->istrip;
169			};
170			if ( plane > splane[view]+gapth ){
171			wgap[view] = true;
172			interpl[view] = splane[view];
173			interv[view] = sview[view];
174			// ii = L2->GetCaloLevel1()->istrip;
175			};
176			splane[view] = plane;
177			sview[view] = view;
178			nhit[view] = 1;
179			} else {
180			nhit[view]++;
181			};
182			};
183			//
184			ii++;
185			//
186			};
187			//
188			// printf("conversion interpl %i interv %i multhit %i interplane %i \n",interpl[0],interv[0],multhit,interplane);
189			Int_t winterplane[2] = {-1,-1};
190			//
191			for ( Int_t view = 0; view < 2; view++){
192			//
193			if ( nhit[view] > 1 && !wmulthit[view] && !wgap[view] ){
194			wmulthit[view] = true;
195			interpl[view] = splane[view];
196			interv[view] = sview[view];
197			};
198			//
199			if ( wmulthit[view] ) multhit = true;
200			if ( wgap[view] ) gap = true;
201			//
202			// convert view and plane number of interaction plane into number of available dE/dx measurements before the interaction plane
203			//
204			if ( interpl[view] >= 0 ) {
205			if ( interv[view] == 0 ){
206			winterplane[view] = (1 + interpl[view]) * 2;
207			} else {
208			winterplane[view] = (1 + interpl[view]) + (1 + interpl[view] - 1);
209			};
210			if ( wmulthit[view] ) winterplane[view]--;
211			};
212			};
213			if ( winterplane[0] > 0 && winterplane[1] > 0 ){
214			if ( multhit ){
215			interplane = TMath::Min(winterplane[0],winterplane[1]);
216			} else {
217			interplane = TMath::Max(winterplane[0],winterplane[1]);
218			};
219			} else {
220			if ( !winterplane[0] \|\| !winterplane[1] ){
221			interplane = 0;
222			} else {
223			interplane = TMath::Max(winterplane[0],winterplane[1]);
224			};
225			};
226			//
227			// printf("2conversion interpl %i interv %i multhit %i interplane %i \n",interpl[1],interv[1],multhit,interplane);
228			// printf("3conversion winterpl0 %i winterpl1 %i \n",winterplane[0],winterplane[1]);
229			//
230			Int_t ipl = 0;
231			if ( interplane > 0 ){
232			//
233			// Calculate preq, postq, qpremean
234			//
235			ii = 0;
236			Int_t ind = -1;
237			Int_t qsplane = -1;
238			Int_t qsview = -1;
239			Float_t qme[200];
240			memset(qme,0,200*sizeof(Float_t));
241			//
242			while ( ii<L2->GetCaloLevel1()->istrip ){
243			//
244			mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);
245			//
246			if ( strip != -1 ){
247			if ( view == 0 ){
248			ipl = (1 + plane) * 2;
249			} else {
250			ipl = (1 + plane) + (1 + plane - 1 );
251			};
252			if ( ipl > interplane ){
253			postq += mip;
254			} else {
255			preq += mip;
256			if ( strip == (track->GetCaloTrack()->tibar[plane][view]-1) \|\| strip == (track->GetCaloTrack()->tibar[plane][view]-2) \|\| strip == (track->GetCaloTrack()->tibar[plane][view]) ){
257			if ( qsplane != plane \|\| qsview != view ){
258			qsplane = plane;
259			qsview = view;
260			ind++;
261			if ( ind > 199 ) printf(" AAAGH!! \n");
262			qme[ind] = 0.;
263			};
264			qme[ind] += mip;
265			};
266			};
267			//
268			};
269			//
270			ii++;
271			//
272			};
273			//
274			// here we must calculate qpremean, order vector qme, select 3 lowest measurements and caculate the mean...
275			//
276			// for (Int_t l=0; l < interplane; l++){
277			// printf(" qme[%i] = %f \n",l,qme[l]);
278			// // if ( qme[ind] < 6.25 ) qme[ind] = 10000.;
279			// };
280			Long64_t work[200];
281			ind = 0;
282			Int_t l = 0;
283			Float_t qm = 0.;
284			Int_t uplim = TMath::Max(3,N);
285			while ( l < uplim && ind < interplane ){
286			qm = TMath::KOrdStat(interplane,qme,ind,work);
287			if ( qm >= mesethr*0.9 ){
288			if ( l < 3 ) qpremean += qm;
289			if ( l < N ) qpremeanN += qm;
290			l++;
291			// printf(" value no %i qm %f \n",l,qm);
292			};
293			ind++;
294			};
295			//
296			qpremean /= l;
297			qpremeanN /= N;
298			//
299			// printf(" charge is %f \n",sqrt(qpremean));
300			//
301			// printf("preq postq\n");
302			//
303			if ( mesethr != ethr && interplane >= 3 && !aldone ){
304			Float_t mesethr2 = (sqrt(qpremean) - 0.50)*(sqrt(qpremean) - 0.50);
305			if ( mesethr2 < mesethr*0.90 ){
306			mesethr = (sqrt(dedx1) - 0.25)*(sqrt(dedx1) - 0.25);
307			} else {
308			mesethr = mesethr2;
309			};
310			aldone = true;
311			if ( mesethr > defethr ) goto retry;
312			};
313			};
314			//
315			// printf(" esci \n");
316			//
317			};