--- calo/flight/CaloNuclei/src/CaloNuclei.cpp	2007/04/03 15:55:16	1.1.1.1
+++ calo/flight/CaloNuclei/src/CaloNuclei.cpp	2008/03/12 14:19:21	1.12
@@ -20,14 +20,23 @@
   PKT = 0;
   atime = 0;
   N = 5;
+  R = 3;
+  //
+  debug = false;
+  usetrack = true;
   //
 };
 
 void CaloNuclei::Clear(){
   //
+  UN = 0;
+  tr = 0;
+  sntr = 0;
   interplane = 0;
   preq = 0.;
   postq = 0.;
+  stdedx1 = 0.;
+  ethr = 0.;
   dedx1 = 0.;
   dedx3 = 0.;
   qpremean = 0.; 
@@ -40,20 +49,24 @@
 
 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 Use track %i \n",OBT,PKT,atime,tr,usetrack);
+  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(" stdedx1 [dE/dx from the first calorimeter plane standalone]: %f \n",stdedx1);
+  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 +77,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");
@@ -74,11 +91,12 @@
   Bool_t newentry = false;
   //
   if ( L2->IsORB() ){
-    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 ){
       newentry = true;
       OBT = L2->GetOrbitalInfo()->OBT;
       PKT = L2->GetOrbitalInfo()->pkt_num;
       atime = L2->GetOrbitalInfo()->absTime;
+      sntr = ntr;
     };
   } else {
     newentry = true;
@@ -86,17 +104,83 @@
   //
   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);
+  if ( debug ) printf(" Always calculate stdedx1 \n");
+  //
+  // Always calculate stdedx1
   //
   Int_t view = 0;
   Int_t plane = 0;
   Int_t strip = 0;
+  Int_t indx = 0;
+  Float_t vfpl[96];
+  Int_t stfpl[96];
+  memset(vfpl, 0, 96*sizeof(Float_t));
+  memset(stfpl, 0, 96*sizeof(Int_t));
   Float_t mip = 0.;
+  for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
+    //
+    mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
+    //
+    // put in vfpl vector the energy release on the first plane
+    //
+    if ( strip != -1 && view == 1 && plane == 0 ) {
+      stfpl[indx] = strip;
+      vfpl[indx] = mip;
+      indx++;
+    };
+    //
+  };
+  //
+  if ( debug ) printf(" find energy released along the strip of maximum on the first plane and on the two neighbour strips \n");
+  //
+  // find energy released along the strip of maximum on the first plane and on the two neighbour strips
+  //
+  if ( indx > 0 ){
+    Int_t mindx = (Int_t)TMath::LocMax(indx,vfpl);
+    for (Int_t ii=0; ii<indx; ii++){
+      if ( stfpl[ii] == stfpl[mindx] ) stdedx1 += vfpl[ii];
+      if ( (mindx-1)>=0 && stfpl[ii] == stfpl[mindx-1] ) stdedx1 += vfpl[ii];
+      if ( (mindx+1)<96 && stfpl[ii] == stfpl[mindx+1] ) stdedx1 += vfpl[ii];
+    };
+  } else {
+    stdedx1 = 0.;
+  };
+  //
+  if ( debug ) printf(" if ( !usetrack ) return: usetrack %i ntr %i \n",usetrack,ntr);
+  // 
+  //
+  //  if ( !usetrack ) return;
+  //
+  PamTrack *ptrack = 0;
+  CaloTrkVar *track = 0;
+  //
+  if ( usetrack ){
+    if ( ntr >= 0 ){
+      ptrack = L2->GetTrack(ntr);
+      if ( ptrack ) track = ptrack->GetCaloTrack();
+    } else {
+      track = L2->GetCaloStoredTrack(ntr);
+    };
+    //
+    if ( !track && ntr >= 0 ){
+      printf(" ERROR: cannot find any track!\n");
+      printf(" ERROR: CaloNuclei variables not completely filled \n");
+      return;   
+    };
+  } else {
+    if ( ntr >= 0 ){
+      if ( debug ) printf(" ERROR: you asked not to use a track but you are looking for track number %i !\n",ntr);
+      if ( debug ) printf(" ERROR: CaloNuclei variables not completely filled \n");
+      return;       
+    };
+  };
+  //
   //  Float_t defethr = 6. * 0.90;
   Float_t defethr = 6.25; // = (sqrt(9) - 0.5) ** 2.;
   //
@@ -106,21 +190,42 @@
     //
     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;
+    if ( ntr >= 0 ){
+      //
+      if ( strip != -1 && 
+	   view == 1   && 
+	   plane == 0  &&
+	   ( strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) || strip == (track->tibar[0][1]) ) 
+	   && true ){      
+	dedx1 += mip;
+      };
+      if ( strip != -1 && 
+	   (( view == 1 && ( plane == 0 || plane == 1 ) ) ||
+	    ( view == 0 && plane == 0 )) &&
+	   (( view == 0 && ( strip == track->tibar[0][0] || strip == (track->tibar[0][0]-1) || strip == (track->tibar[0][0]-2) )) ||
+	    ( view == 1 && ( strip == track->tibar[0][1] || strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) )) ||
+	    ( view == 1 && ( strip == track->tibar[1][1] || strip == (track->tibar[1][1]-1) || strip == (track->tibar[1][1]-2) ))) &&
+	   true ){
+	dedx3 += mip;
+      };
+    } else {
+      //
+      if ( strip != -1 && 
+	   view == 1   && 
+	   plane == 0  &&
+	   ( strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) || strip == (L2->GetCaloLevel2()->cibar[0][1]) ) 
+	   && true ){      
+	dedx1 += mip;
+      };
+      if ( strip != -1 && 
+	   (( view == 1 && ( plane == 0 || plane == 1 ) ) ||
+	    ( view == 0 && plane == 0 )) &&
+	   (( view == 0 && ( strip == L2->GetCaloLevel2()->cibar[0][0] || strip == (L2->GetCaloLevel2()->cibar[0][0]-1) || strip == (L2->GetCaloLevel2()->cibar[0][0]-2) )) ||
+	    ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[0][1] || strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) )) ||
+	    ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[1][1] || strip == (L2->GetCaloLevel2()->cibar[1][1]-1) || strip == (L2->GetCaloLevel2()->cibar[1][1]-2) ))) &&
+	   true ){
+	dedx3 += mip;
+      };
     };
     //    
   };
@@ -133,7 +238,7 @@
   //
  retry:
   // 
-  //  printf("retry\n");
+  if ( debug ) printf("retry\n");
   //
   interplane = 0;
   //
@@ -154,30 +259,63 @@
     //
     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 ( ntr >= 0 ){
+      if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
+	   ( strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) ) 
+	   && true ){      
+	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]);
+	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];	  
 	};
-	if ( plane > splane[view]+gapth ){
-	  wgap[view] = true;
-	  interpl[view] = splane[view];
-	  interv[view] = sview[view];
-	  //	  ii = L2->GetCaloLevel1()->istrip;
+	  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];
+	  };
+	  splane[view] = plane;
+	  sview[view] = view;
+	  nhit[view] = 1;
+	} else {
+	  nhit[view]++;
+	};
+      };
+    } else {
+      if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
+	   ( strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) ) 
+	   && true ){      
+	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]);
+	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];	  
+	};
+	  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];
+	  };
+	  splane[view] = plane;
+	  sview[view] = view;
+	  nhit[view] = 1;
+	} else {
+	  nhit[view]++;
 	};
-	splane[view] = plane;
-	sview[view] = view;
-	nhit[view] = 1;
-      } else {
-	nhit[view]++;
       };
     };
     //
@@ -185,7 +323,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 +362,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 +374,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 ){
       //
@@ -253,15 +396,54 @@
 	  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.;
+	  if ( ntr >= 0 ){
+	    if (  strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) ){
+	      if ( qsplane != plane || qsview != view ){
+		qsplane = plane;
+		qsview = view;
+		ind++;
+		if ( debug && ind > 199 ) printf(" AAAGH!! \n");
+		qme[ind] = 0.;
+	      };
+	      qme[ind] += mip;
 	    };
-	    qme[ind] += mip;
+	    for ( Int_t ns = 0; ns < R ; ns++){
+	      Int_t ms =  track->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;
+	      };
+	    };	
+	  } else {
+	    if (  strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) ){
+	      if ( qsplane != plane || qsview != view ){
+		qsplane = plane;
+		qsview = view;
+		ind++;
+		if ( debug && ind > 199 ) printf(" AAAGH!! \n");
+		qme[ind] = 0.;
+	      };
+	      qme[ind] += mip;
+	    };
+	    for ( Int_t ns = 0; ns < R ; ns++){
+	      Int_t ms =  L2->GetCaloLevel2()->cibar[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 +455,58 @@
     //
     // 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;
+    Int_t RN = 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 ( l < 3 ) qpremean += qm;
-	if ( l < N ) qpremeanN += qm;
+      if ( qm >= qmt ){	
+	if ( l < 3 ){
+	  qpremean += qm;
+	  RN++;
+	};
 	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;
+    qpremean /= (Float_t)RN;
     //
-    //  printf(" charge is %f \n",sqrt(qpremean));
+    ind = 0;
+    l = 0;
+    RN = 0;
+    while ( l < uplim && ind < interplane ){
+      qm2 = TMath::KOrdStat(interplane,qme2,ind,work);
+      if ( qm2 >= qmt ){	
+	if ( l < N ){
+	  qpremeanN += qm2;
+	  RN++;
+	};
+	l++;
+	if ( debug ) printf(" qm2 value no %i qm %f qmt %f RN %i \n",l,qm2,qmt,RN); 
+      };
+      ind++;
+    };
     //
-    //    printf("preq postq\n");    
+    qpremeanN /= (Float_t)RN;
+    UN = RN;
+    //
+    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);
@@ -308,10 +516,20 @@
 	mesethr = mesethr2;
       };
       aldone = true;
-      if ( mesethr > defethr ) goto retry;
+      if ( mesethr > defethr ){
+	interplane = 0;
+	preq = 0.;
+	postq = 0.;
+	qpremean = 0.; 
+	qpremeanN = 0.; 
+	multhit = false;
+	gap = false;
+	goto retry;
+      };
     };
   };
   //
-  //  printf(" esci \n");
+  if ( debug ) this->Print();
+  if ( debug ) printf(" esci \n");
   //
 };