--- calo/flight/CaloNuclei/src/CaloNuclei.cpp	2007/04/03 15:55:16	1.1
+++ calo/flight/CaloNuclei/src/CaloNuclei.cpp	2007/04/05 13:38:32	1.4
@@ -20,11 +20,15 @@
   PKT = 0;
   atime = 0;
   N = 5;
+  R = 3;
+  //
+  debug = false;
   //
 };
 
 void CaloNuclei::Clear(){
   //
+  tr = 0;
   interplane = 0;
   preq = 0.;
   postq = 0.;
@@ -40,20 +44,23 @@
 
 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");
+  Process();
+  //
+  printf("========================================================================\n");
+  printf(" OBT: %u PKT: %u ATIME: %u Track %i \n",OBT,PKT,atime,tr);
+  printf(" interplane [number of available dE/dx before interaction]: %i\n",interplane);
+  printf(" ethr [threshold used to determine interplane]:............ %f \n",ethr); 
+  printf(" dedx1 [dE/dx from the first calorimeter plane]:........... %f \n",dedx1);
+  printf(" dedx3 [dE/dx (average) if the first 3 Si planes]:......... %f \n",dedx3);
+  printf(" multhit [true if interplane determined by multiple hits]:. %i \n",multhit);
+  printf(" gap [true if interplane determined by a gap]:............. %i \n",gap);
+  printf(" preq [total energy in MIP before the interaction plane]:.. %f \n",preq);
+  printf(" postq [total energy in MIP after the interaction plane]:.. %f \n",postq);
+  printf(" qpremean [truncated mean using 3 planes and 3 strips]:.... %f \n",qpremean); 
+  printf(" N [no of used plane]:..................................... %i \n",N); 
+  printf(" R [no strip used per plane ]:............................. %i \n",R); 
+  printf(" qpremeanN [truncated mean using N planes and R strips]:... %f \n",qpremeanN); 
+  printf("========================================================================\n");
   //
 };
 
@@ -64,7 +71,11 @@
 
 
 void CaloNuclei::Process(){
-  //
+  Process(0);
+};
+
+void CaloNuclei::Process(Int_t ntr){
+  //  
   if ( !L2 ){
     printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
     printf(" ERROR: CaloNuclei variables not filled \n");
@@ -86,12 +97,14 @@
   //
   if ( !newentry ) return;
   //
-  //  printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
+  tr = ntr;
+  //
+  if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
   //
   Clear();
   //
   PamTrack *track = 0;
-  track = L2->GetTrack(0);
+  track = L2->GetTrack(ntr);
   //
   Int_t view = 0;
   Int_t plane = 0;
@@ -133,7 +146,7 @@
   //
  retry:
   // 
-  //  printf("retry\n");
+  if ( debug ) printf("retry\n");
   //
   interplane = 0;
   //
@@ -157,21 +170,23 @@
     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]);
+      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]);
       interpl[view] = plane;
       interv[view] = view;
       if ( splane[view] != plane || sview[view] != view ){
 	if ( nhit[view] > 1 ){
 	  wmulthit[view] = true;
+	  //	  if ( splane[view] == -1 ) splane[view] = 0; //
+	  //	  if ( sview[view] == -1 ) sview[view] = view; //
 	  interpl[view] = splane[view];
-	  interv[view] = sview[view];
-	  // 	  ii = L2->GetCaloLevel1()->istrip;
+	  interv[view] = sview[view];	  
 	};
 	if ( plane > splane[view]+gapth ){
 	  wgap[view] = true;
+	  //	  if ( splane[view] == -1 ) splane[view] = 0;//
+	  //	  if ( sview[view] == -1 ) sview[view] = view; //
 	  interpl[view] = splane[view];
 	  interv[view] = sview[view];
-	  //	  ii = L2->GetCaloLevel1()->istrip;
 	};
 	splane[view] = plane;
 	sview[view] = view;
@@ -185,7 +200,7 @@
     //    
   };
   //
-  //  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);
   Int_t winterplane[2] = {-1,-1};
   //
   for ( Int_t view = 0; view < 2; view++){
@@ -224,8 +239,8 @@
     };
   };
   //  
-  //  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]);
+  if ( debug ) printf("2conversion interpl %i interv %i multhit %i interplane %i \n",interpl[1],interv[1],multhit,interplane);
+  if ( debug ) printf("3conversion winterpl0 %i winterpl1 %i \n",winterplane[0],winterplane[1]);
   //
   Int_t ipl = 0;
   if ( interplane > 0 ){
@@ -236,8 +251,13 @@
     Int_t ind = -1;
     Int_t qsplane = -1;
     Int_t qsview = -1;
+    Int_t ind2 = -1;
+    Int_t qsplane2 = -1;
+    Int_t qsview2 = -1;
     Float_t qme[200];
     memset(qme,0,200*sizeof(Float_t));
+    Float_t qme2[2112];
+    memset(qme2,0,2112*sizeof(Float_t));
     //
     while ( ii<L2->GetCaloLevel1()->istrip ){
       //
@@ -258,11 +278,24 @@
 	      qsplane = plane;
 	      qsview = view;
 	      ind++;
-	      if ( ind > 199 ) printf(" AAAGH!! \n");
+	      if ( debug && ind > 199 ) printf(" AAAGH!! \n");
 	      qme[ind] = 0.;
 	    };
 	    qme[ind] += mip;
 	  };
+	  for ( Int_t ns = 0; ns < R ; ns++){
+	    Int_t ms =  track->GetCaloTrack()->tibar[plane][view] - 1 - ns + (R - 1)/2;
+	    if ( strip == ms ){
+	      if ( qsplane2 != plane || qsview2 != view ){
+		qsplane2 = plane;
+		qsview2 = view;
+		ind2++;
+		if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
+		qme2[ind2] = 0.;
+	      };
+	      qme2[ind2] += mip;
+	    };
+	  };	
 	};      
 	//
       };
@@ -273,32 +306,49 @@
     //
     // 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.;      
-    //     };
+    if ( debug ){
+      for (Int_t l=0; l < interplane; l++){
+	printf(" qme[%i] = %f \n",l,qme[l]);  
+      };
+    };
+    //
     Long64_t work[200];
     ind = 0;
     Int_t l = 0;
     Float_t qm = 0.;
+    Float_t qm2 = 0.;
+    //
+    Float_t qmt = ethr*0.8; // *0.9
+    //
     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 ( qm >= qmt ){	
 	if ( l < 3 ) qpremean += qm;
-	if ( l < N ) qpremeanN += qm;
 	l++;
-	//	printf(" value no %i qm %f \n",l,qm); 
+	if ( debug ) printf(" value no %i qm %f qmt %f \n",l,qm,qmt); 
       };
       ind++;
     };
     //
     qpremean /= l;
-    qpremeanN /= N;
     //
-    //  printf(" charge is %f \n",sqrt(qpremean));
+    ind = 0;
+    l = 0;
+    while ( l < uplim && ind < interplane ){
+      qm2 = TMath::KOrdStat(interplane,qme2,ind,work);
+      if ( qm2 >= qmt ){	
+	if ( l < N ) qpremeanN += qm2;
+	l++;
+	if ( debug ) printf(" qm2 value no %i qm %f qmt %f \n",l,qm2,qmt); 
+      };
+      ind++;
+    };
+    //
+    qpremeanN /= l;
     //
-    //    printf("preq postq\n");    
+    if ( debug ) printf(" charge is %f \n",sqrt(qpremean));
     //
     if ( mesethr != ethr && interplane >= 3 && !aldone ){
       Float_t mesethr2 = (sqrt(qpremean) - 0.50)*(sqrt(qpremean) - 0.50);
@@ -312,6 +362,6 @@
     };
   };
   //
-  //  printf(" esci \n");
+  if ( debug ) printf(" esci \n");
   //
 };