// This class handels the voltage scans produced with TCT Labview software (".tct" files)!
// Waveform vectors are stored in the histograms on which the data manipultaion rutines work
#include "Rtypes.h"
#include <math.h>
//#include <Riostream>
#include <stdio.h>
#include <stdlib.h>
#include "Analog.h"
#include <string.h>
#include "TF1.h"
#include "TDatime.h"
#include "TCanvas.h"

#define MXREG  50

#ifdef TIMING
float timex_(float *);
float p01,p02,p11,p22,p1,p2,ptim[10],xtim[10];
#define INITIMING { float t0; timest_(&t0); }
#define STTIME(p)   { timex_(&p); }
#define EDTIME(i,p1,p2) { 			\
  timex_(&p2); 					\
  ptim[i]+=p2-p1;				\
  xtim[i]++;                                    \
}
#else
#define INITIMING
#define STTIME(p)
#define EDTIME(i,p1,p2)
#endif

#define  SCALE_NOISE    0.89
#define  NSNAP          10
#define  NCHAN_DEFAULT 128
#define  my_way1 0 /* Calculate average noise within the program */
#define  my_way2 0 /* fill histo signal/aver_noise, yes=1, no=0 */
#define  my_way3 0 /* scale noise yes=1, no=0 */
#define  SUBcmod(x,i,j) (x - cmod[j][0])



void swoo(char *a, char *b) {
  char c = *a;
  *a = *b; *b = c;
}

void swooip(float *in, int s) {
  char *sr, b;
  while(s--) {
    sr=(char *)in;
    swoo(&sr[0], &sr[3]);
    swoo(&sr[1], &sr[2]);
    in++;
  }
}               


ClassImp(Analog)


//____________________
// MeassureWF
//                                                             //
// This class handels the voltage scans produced with TCT Labview software (".tct" files)!
// Waveform vectors are stored in the histograms on which the data manipultaion rutines work

//Analog::Analog(Char_t *name,Int_t num, Int_t start,Int_t bins,Float_t startbin,Float_t endbin)

 Analog::Analog(char *filename,int numch)
{

   register int i,j;
   int      nevt;
   
   NNchan=numch;
   lrecl = NNchan*sizeof(float);
   chan  = new CHannel [NNchan];  
   for(i=0;i<NNchan;i++)  {chan[i].id=i+1;}

   for(i=0;i<200;i++) his[i]=NULL;
   for(i=0;i<5;i++)   his2[i]=NULL;
   for(i=0;i<10;i++)  {prof[i]=NULL; chis[i]=NULL;}


/*------------------------------------------------------------------------
 * Open data file
 *------------------------------------------------------------------------*/

   strcpy(ffile,filename);
   if((in=fopen(ffile,"rb"))==NULL) 
     {
      printf("Failed to open the data filen"); exit(0);
     }
   
   fseek(in,0,SEEK_END);
   Nevt = ftell(in)/lrecl;
   fseek(in,0,SEEK_SET);

   sample = (sample >0 ? Nevt/sample : 0 );
   printf("=> About %d events in the filen",Nevt);
   
}


 void Analog::pedes(int npass,int se,int ee)
{
/*@(#)******************************************************************
 * pedes
 *
 *   Purpose: Pedestal calculation
 *   Inputs :
 *        npass - iteration number passing through data (=1 no threshold)
 *        nchan - number of channels
 *        nregs - no. of defined regions
 *        regs  - pointer to defined regions
 *        data  - data
 *
 **********************************************************************/
   register int i,ii,j,jj;
   register CHannel *pdat;
   int      ievt,*noise_dead;
   double   *dbuff,*sum,*sumsq,*xn,ddum,sdum;
   CHannel  *pend;
   Region   *preg;
   float    *buff,*ptr,fdum,fsgn,cut;
   double    *avnoidead;
   CM       *cmod;

/*------------------------------------------------------------------------
 * Get space to compute mean value and deviation
 *------------------------------------------------------------------------*/
   noise_dead=(int *)malloc(sizeof(int)*NNchan);
   avnoidead=(double *)malloc(sizeof(double)*Nreg);
   avnoise=(double *) calloc(Nreg,sizeof(double));
   dbuff = (double *) calloc(3*NNchan,sizeof(double));

   sum   = dbuff;
   sumsq = sum   + NNchan;
   xn    = sumsq + NNchan;
   buff = (float *)  calloc(NNchan,sizeof(float));
   cmod = (CM *) calloc(Nreg,sizeof(CM));
   if (sample<=0) sample=1;
/*------------------------------------------------------------------------
 * Rewinds the file and starts reading
 *------------------------------------------------------------------------*/
   //   prepare_file(pfile,0,SEEK_SET);
   //prepare_file(in,0,SEEK_SET);

   prepare_file(in,lrecl*se,SEEK_SET);

   pend = chan+NNchan;
   ievt = 0; 
   //   while (GetEvent(pfile,buff,lrecl)) {
   while (GetEvent(buff,lrecl) && (ee>0?ievt<ee:true)  ) {
     ievt++; //printf("npass=%d event %dn",npass,se+ievt);
      /*  Get data   */
      for (pdat=chan,ptr=buff;pdat<pend;pdat++,ptr++) {
	//	printf("To je stevilka,%f,n",*ptr); 

	 pdat->raw = ( pdat->excl ? 999999. : (1.+ *ptr)*1000.*ssign);
	 pdat->signal = pdat->raw - pdat->pd;
      }
      /*  Compute Common Mode */
      if (npass) CommonMode(cmod);

      /*  For pedestals calculations */
      for (j=0;j<Nreg;j++) {
	 preg = reg+j;
	 pdat = chan+preg->ifi;
	 for (i=preg->ifi;i<=preg->ila;i++,pdat++) {
	    cut  = ( npass ? COMhard*pdat->sd : 9999.);
	    fdum = SUBcmod(pdat->raw,i,j);
	    fsgn = fdum - pdat->pd;
	    //if (npass==1)  HF1(3800,fsgn,1.);
            if ( fabs(fsgn)<cut) {
	      //if (npass==2) HF1(3700+i,fsgn,1.);
               sum[i]   += fdum;
	       sumsq[i] += fdum*fdum;
	       xn[i]++;
	    }
	 }
      }
      
      if ( sample>1) prepare_file(in,sample*lrecl,SEEK_CUR);
	//prepare_file(pfile,smp*lrecl,SEEK_CUR);
   }

/*------------------------------------------------------------------------
 * Compute pedestals and noise
 *------------------------------------------------------------------------*/

   for (j=0;j<Nreg;j++) {
        
      ii=0; jj=0; avnoidead[j]=0.0; avnoise[j]=0.0;
      preg = reg+j;
      pdat = chan+preg->ifi;
      for (i=preg->ifi;i<=preg->ila;i++,pdat++) {

	 if ( xn[i]> 1 ) {
	    ddum    = sum[i]/xn[i];
	    sdum    = (sumsq[i]-xn[i]*ddum*ddum)/(xn[i]-1.);
	    pdat->pd0 = pdat->pd = (float) ddum;
	    pdat->sd0 = pdat->sd = (float) sqrt(fabs(sdum)); //absolute value is needed because of wrong sqrt(x)
	 }
	 else
	    pdat->pd = pdat->sd = 0;

	 pdat->signal = SUBcmod(pdat->raw,i,j)-pdat->pd;
	 pdat->sn     = (pdat->sd !=0 ? pdat->signal/pdat->sd : 0);

	 //	 printf("%f %fn",pdat->pd,pdat->sd);
	 his[npass*2]->Fill(i+0.5,pdat->pd);
	 his[npass*2+1]->Fill(i+0.5,pdat->sd);
	 //	 if(my_way3 && npass==2) {pdat->sd=pdat->sd*SCALE_NOISE;
	 //	                          HFF1(3850,nid[npass],i+0.5,pdat->sd);}

	 if(pdat->sd > 0.3) {avnoise[j]+=pdat->sd; ii++;}
	                  else
			   {avnoidead[j]+=pdat->sd; jj++; noise_dead[i]=1;
			   if(my_way1) pdat->excl=1;
			   }
      }
    avnoise[j]=avnoise[j]/ii; avnoidead[j]=avnoidead[j]/jj;
    printf("n Reg. #: %d  Avn: %d %f  Dn.: %d %f",j,ii,avnoise[j],jj,avnoidead[j]);
  }
   printf("n");
   free(dbuff);
   free(buff);
   free(cmod);

}


 int Analog::GetEvent(void *buff, unsigned int lrec, int en)
{
 
     
   size_t iok;
   int nc,i,jj,yes;
   fpos_t pos;
   float *pw,*pr,*pend,*pb,*yyy;
   void swooip(float *,int);

  if(en!=-1 && en>0)
     prepare_file(in,lrec*en,SEEK_SET);
   STTIME(p1);
   //   iok = _read(pfile,buff,lrec);
   iok=fread(buff,sizeof(float),lrec/sizeof(float),in); 
   //      printf("iok=%d , %d, gregorn",lrec,iok);
      //yes=fgetpos(in,&pos );
      //if(pos%lrec!=0) printf("Sranje *******************************");
   yyy=(float *)buff;
   if (iok && skip ) {
      pb = (float *) buff;
      pend = pb+NNchan;
      if (skip==1) {
	 pr = pb+2;
	 pw = pb+1;
      }
      else {
	 pr = pb+1;
	 pw = pb;
      }
	

      for ( ;pr<pend;pr+=2,pw++) *pw = *pr;
      for ( ;pw<pend;pw++) *pw=0;
   }
   EDTIME(0,p1,p2);
   swooip(yyy,NNchan);
   //printf("n raw:data Nchan=%d",NNchan);  
   //while(i++<31) if(yyy[i]>10){ 
   //if(pos%(20000)==0)   {for(jj=0;jj<50;jj++) printf(" %e",*(yyy+jj)); printf("  pos=%dn",pos);}
   //   if(pos==2999*128) return(0);
   return (iok);
}


 void Analog::CommonMode(CM *cmod)
{
   register int i,j,ip;
   register CHannel *ptr;
   float sy,sy2,xn,cut,det;
   float  dev;
   Region  *preg;
   float   *cmean,*csig,ccut;

   static   int sCM=sizeof(CM);

   memset(cmod,0.,Nreg*sCM);
   for ( ip=0;ip<2;ip++) {
      for (i=0;i<Nreg;i++) {
	 preg=reg+i;
	 cmean= &cmod[i][0];
	 csig = &cmod[i][1];
	 ccut = (*csig)*COMthres;
	 xn = sy = sy2 = 0;
	 ptr = chan+preg->ifi;
	 for (j=preg->ifi;j<=preg->ila;j++,ptr++) {
	    if ( ptr->excl ) continue;
	    cut = COMhard*ptr->sd;
	    dev = ( ip==0 ? -1. : fabs(ptr->signal- *cmean) );
	    if ( fabs(ptr->signal)<cut && dev<ccut ) {
	       sy  += ptr->signal;
	       sy2 += ptr->signal*ptr->signal;
	       xn++;
	    }	
	 }
	 if (xn>0.) *cmean = sy/xn;
	 if (xn>2.) *csig  = sqrt( ((double)sy2-xn*(*cmean)*(*cmean))/(xn-1.) );
      }
   }
}






 Analog::~Analog()
{
  delete [] chan;
  fclose(in);
}

////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////
 void Analog::SetRegions(int numr,int regions[])
{
  Nreg=numr;
  reg=(Region *) calloc(Nreg,sizeof(Region));
  for(int i=0;i<Nreg;i++)
    {
      reg[i].ifi=regions[i*2];
      reg[i].ila=regions[i*2+1];
    }
}

 void Analog::SetExclude(int nume,int *exclude)
{
  register int i;
  if(nume==0) Nexcl=0; else Nexcl=nume;
  excld = (int *) calloc(Nexcl,sizeof(int));
  for( i=0;i<nume;i++) excld[i]=exclude[i];
  for( i=0;i<Nexcl;i++) chan[excld[i]-1].excl=1; 
}


 void Analog::GetParams()
{ 
  int k,i;
/*------------------------------------------------------------------------
 * Write what we got
 *------------------------------------------------------------------------*/
   printf("n=========================================================n");
   printf("%s - file n",ffile);
   printf("no. of channels/record : %9dn",NNchan);
   printf("Pedestal update rate   : %9dn",nupdate);
   printf("Expected sign of signal: %9.4fn",ssign);
   printf("Hard threshold         : %9.4fn",COMhard);
   printf("Common mode Threshold  : %9.4fn",COMthres);
   printf("Clustering thresholds  :n");
   printf("tSeed     : %9.4fn\tInclusion: %9.4fn\tCluster  : %9.4fn",
	  CLUSthr[0],CLUSthr[1],CLUSthr[2]);
   printf("tmax. no. of clusters & strips /cluster : %d %dn",mxclust,mxstrip);
   printf("Regions analyzed       : %4dn",Nreg);
   for (i=0;i<Nreg;i++) {
      printf("t%3d -> %3dn",reg[i].ifi,reg[i].ila);
      reg[i].ifi--;
      reg[i].ila--;
   }
   printf("Channels Excluded      : %4dn",Nexcl);
   for (i=0,k=0;i<Nexcl;i+=10) {
      printf("t");
      while (k<10) {
	 if ( i+k<Nexcl )
	    printf("%3d ",excld[i+k]);
	 k++;
      }
      printf("n");
      k=0;
   }
   printf("=========================================================n\n");

}



 void Analog::ana(int se,int ee,int histos)
{
   register int i,j,k,l,m;
   int      ievt=0,sevt=0,iflg;
   float    *buff,*ptr;
   CHannel  *pdat,*pend;
   CM       *cmod;
   int      *nclst,mxstrip;
   lCLUster *clst,pcls,pc;
   Region  *preg;
   float    sn,sn0;
   float    par[4],sig[4],chi2;
   char     fnam[4][8]={"Scale","Mean","Norma","Sigma"};

   STTIME(p1);
   buff  = (float *)  calloc(NNchan,sizeof(float));
   cmod  = (CM *) calloc(Nreg,sizeof(CM));
   clst  = (lCLUster *) calloc(Nreg,sizeof(lCLUster));
   nclst = (int *) calloc(Nreg,sizeof(int));
   for (i=0;i<Nreg;i++)
      clst[i] = (CLUster *) calloc(reg->ila - reg->ifi + 1, sizeof(CLUster));
   EDTIME(1,p1,p2);
/*------------------------------------------------------------------------
 * Rewinds the file and starts reading
 *------------------------------------------------------------------------*/
   //   prepare_file(pfile,0,SEEK_SET);
     prepare_file(in,lrecl*se,SEEK_SET);
     //prepare_file(in,0,SEEK_SET);
   pend = chan+NNchan;
   //   while (GetEvent(pfile,buff,lrecl)) {

   if(histos==2) Reset_histograms(1);


   while (GetEvent(buff,lrecl) && (ee>0?ievt<ee:true)) {
      STTIME(p01);
      //      printf("event %dn",se+ievt);      
      ievt++; 
      /*  Get data   */
      STTIME(p1);
      pdat=chan;
      ptr=buff;
      for (i=0;i<NNchan;i++,pdat++,ptr++) {
       	 pdat->raw = ( pdat->excl ? 999999. : (1.+ *ptr)*1000.*ssign);
	 pdat->signal = pdat->raw - pdat->pd;
	 pdat->clus=0;
      }
      EDTIME(2,p1,p2);
      STTIME(p1);
      /*  Compute Common Mode */
      CommonMode(cmod);
      EDTIME(3,p1,p2);
      /* Loop on strips */
      STTIME(p11);
      iflg = 0;
      for (j=0;j<Nreg;j++) {
	 preg=reg+j;
         k=j+1;
/*------------------------------------------------------------------------
 * Remove pedestals and cm
 *------------------------------------------------------------------------*/
	 STTIME(p1);

	 if(histos)
	   {
	     his[20*k]->Fill(cmod[j][0],1.);
	     his[20*k+1]->Fill(cmod[j][1],1.);
	   }
	 for (i=preg->ifi;i<=preg->ila;i++) {
	    pdat = chan+i;
	    sn   = SUBcmod(pdat->raw,i,j);
	    sn0  = sn - pdat->pd;
	    /*  Update pedestals */
	    if (nupdate > 0. && ievt> (int) nupdate && fabs(sn0)<COMhard*pdat->sd) {
	       pdat->pd = (pdat->pd + sn/nupdate)/(1+1./nupdate);
	       sn -= pdat->pd;
	       pdat->sd = (float) sqrt((pdat->sd*pdat->sd + (sn*sn)/nupdate)/(1+1./nupdate)) ;

	   if(histos)
	    {
	       prof[0]->Fill((float) i,pdat->pd,1.);
	       prof[1]->Fill((float) i,pdat->sd,1.);
	    }
	    }
	    pdat->signal = sn - pdat->pd;
	    pdat->sn     = (pdat->sd !=0 ? pdat->signal/pdat->sd : 0);
	     if(histos) his2[0]->Fill((float) i,SUBcmod(pdat->raw,i,j),1.);
	    }
	 EDTIME(4,p1,p2);

	 /*----------------------------------------------------------------
	  * check stabilities
	  *---------------------------------------------------------------*/
	 /*
	 HCDIR("STAB"," ");	
	 for (i=preg->ifi;i<=preg->ila;i++) {
	    pdat = data+i;
	    sn = SUBcmod(pdat->raw,i,j);
	    sn0 = (sn - pdat->pd);
	    s2  = (sn - pdat->pd0);
	    HFILL(i+1,(float) ievt,sn,1.);
	    HFILL(1001+i,(float) ievt,pdat->pd,1.);
	    HFILL(2001+i,(float) ievt,pdat->sd0,1.);
	    HFILL(3001+i,(float) ievt,pdat->sd,1.);
	 }
	 HF1(4000,pdat->sd0-pdat->sd,1.);
	 HCDIR("\"," ");
	 */

/*------------------------------------------------------------------------
 * Search for clusters
 *------------------------------------------------------------------------*/
	 STTIME(p1);
         pcls = clst[j];
	 nclst[j]=Cluster(preg->ila-preg->ifi+1,chan+preg->ifi,
			  CLUSthr[0],CLUSthr[1],CLUSthr[2],
			  pcls,&mxstrip);
	 if ( nclst[j]>1 || mxstrip>=5 ) continue;
/*------------------------------------------------------------------------
 * Fill histograms
 *------------------------------------------------------------------------*/
	 his[20*k+6]->Fill((float) nclst[j],1.);
	 if ( nclst[j]>0 ) {
	    iflg++;
	    for (i=0,pc=pcls;i<nclst[j];i++,pc++) {
	       l = 1+(pc->fst-chan);
	       m =  (pc->hig-chan) +1;
	       his[20*k+2]->Fill(pc->sn,1.);
	       //	       if(my_way2) HF1(3900+k,pc->sig/avnoise[k],1.);
                if(histos)
		  {
		    his[20*k+3]->Fill(pc->sig,1.);
		    his[20*k+4]->Fill(pc->hig->sn,1.);
		    his[20*k+5]->Fill(pc->hig->signal,1.);
		    his[20*k+7]->Fill((float) pc->ns,1.);
		    his[20*k+8]->Fill((float)pc->eta,1.);
		    if (pc->ns==1) { his[20*k+9]->Fill((float) pc->sn,1.);
	                        his[20*k+10]->Fill((float)pc->sig,1.);}  
		    his[7]->Fill((float) m,1.);
		    prof[2]->Fill((float) m,pc->hig->signal,1.);
		  }
	    }
	    } 
	 EDTIME(5,p1,p2);
      }
      EDTIME(6,p11,p22);
   }


/*------------------------------------------------------------------------
 * Fit s/n plots
 *------------------------------------------------------------------------*/
//    printf("Done. Fitting now S/N distributionsn");
//    for (j=1;j<Nreg+1;j++) {
//       i = 3000+j;
//       par[0] = 1.;
//       par[1] = 0; //HSTATI(i,1," ",0);
//       par[2] = 0; //HMIN(i);
//       par[3] = 2.;
//       chi2   = 0.;
//       HFITH(i,landau,"FQU",4,par,NULL,NULL,NULL,sig[0],chi2);
//       printf("tRegion %d : s/n --> %9.4f width ---> %9.4fn",
// 	     i,par[1]+0.225*par[0],par[3]);
//       HFINAM(i,fnam,4);
//    }

   free(buff);
   free(cmod);
   for (i=0;i<Nreg;i++) free(clst[i]);
   free(clst);
   free(nclst);
}




/*@(#)*****************************************************************
 * prepare_file
 *
 *   Created:  17.9.2007        Author: Gregor Kramberger
 *   Purpose:
 *
 **********************************************************************/
 void Analog::prepare_file(FILE *in,off_t offset, int whence)
{
  //_lseek(pfile,offset,whence);
      fseek(in, offset, whence );
}

/*@(#)*****************************************************************
 *   Book_histograms
 *
 *   Created:  17.9.2007        Author: Gregor Kramberger
 *   Purpose:
 *
 **********************************************************************/

 void Analog::Book_histograms(int nevt)
{
   float  *pdat,*pend,*pin,xsig;
   register int i,ii;
   int iok;
   float mnsig,mxsig;

   his[0]=new TH1F("pedestals1","pedestals, 1 pass",128,.5,128.5);
   his[1]=new TH1F("sigmas1","sigmas, 1 pass",128,.5,128.5);
   his[2]=new TH1F("pedestals2","pedestals, 2 pass",128,.5,128.5);
   his[3]=new TH1F("sigmas2","sigmas, 2 pass",128,.5,128.5);
   his[4]=new TH1F("pedestals3","pedestals, 3 pass",128,.5,128.5);
   his[5]=new TH1F("sigmas3","sigmas, 3 pass",128,.5,128.5);
   prof[0]=new TProfile("ped. update","ped.(update)",128,.5,128.5,-9999.,9999.,"S");
   prof[1]=new TProfile("sigmas (updated)","sigmas (updated)",128,.5,128.5,-9999.,9999.,"S");
   his[7]=new TH1F("hit channel","hit channel",128,.5,128.5);

   for (i=1;i<=Nreg; i++) 
     {
      his[20*i]=new TH1F("cmod mean","cmod mean",96,-20.,20.);
      his[20*i+1]=new TH1F("cmod sigma","cmod sigma",50,-5.,5.);
      his[20*i+2]=new TH1F("cluster S/N","cluster S/N",100,0.,50.);
      //      if(my_way2) HBOOK1(3900+i,"signal over aver. noise",100,0.,50.,0.);
      his[20*i+3]=new TH1F("cluster signal (mV)","cluster signal (mV)",100,0.,150.);
      his[20*i+4]=new TH1F("cluster S/N 1","cluster S/N 1",100,0.,50.);
      his[20*i+5]=new TH1F("cluster signal 1 (mV)","cluster signal 1 (mV)",100,0.,150.);
      his[20*i+6]=new TH1F("Number of clusters","Number of clusters",10,0.,10.);
      his[20*i+7]=new TH1F("no. of strips per cluster","no. of strips per cluster",10.,0.,10.);
      his[20*i+8]=new TH1F("eta","eta",50,0.,1.);
      his[20*i+9]=new TH1F("1 strip S/N","1 strip S/N",50,0.,50.);
      his[20*i+10]=new TH1F("1 strip sig","1 strip sig",50,0.,100.);
     }

   pdat  = (float *)  calloc(NNchan,sizeof(float));
   mnsig= 99999.;
   mxsig=-99999.;
   /* Getr maximum and minimum value of signal */
   //   iok = _read(pfile,pdat,lrecl);
   iok=fread(pdat,sizeof(float),lrecl/sizeof(float),in); 
   for (pin=pdat,pend=pdat+iok/sizeof(float);pin<pend;pin++) {
      xsig = (1.+ *pin)*1000.*ssign;
      mnsig = ( mnsig>xsig ? xsig : mnsig);
      mxsig = ( mxsig<xsig ? xsig : mxsig);
   }
   xsig  = (mxsig-mnsig);
   mnsig -= 0.1*xsig;
   mxsig += 0.1*xsig;

   his2[0]=new TH2F("Channel signal","Channel signal",128,.5,128.5,75,mnsig,mxsig);
   //   HBSLIY(100,128,0.);
   prof[2]=new TProfile("Channel average signal","Channel average signal",128,.5,128.5,-9999.,9999.," ");

   free(pdat);
}

/*@(#)******************************************************************
 * Clus
 *
 *   Created: Thu Nov  7 23:22:50 1996        Author: Carlos Lacasta
 *   Purpose: Searches for Clusters
 *   Input:
 *            ns      - number of strips to look at
 *            data    - pointer to first strip in the list
 *            SeedT   - threshold for cluster seed
 *            IncT    - threshold for neighbours
 *            clsT    - threshold for cluster
 *   Output:
 *            mxstrip - mx. no. of strips found in a cluster
 *            clst    - pointer to an array of clusters
 *   Returns:
 *            number of clusters found
 *
 **********************************************************************/
 int Analog::Cluster(int ns, CHannel *data, float SeedT, float IncT, float clsT,
	 CLUster *clst, int *mxstrip)
{
   int nstrip,nleft,nright,nclus;
   CHannel *ptr=data,*pend=data+ns,*pr,*pl,*ph;
   CLUster *pclst=clst;
   float   sn,sg;

   nclus=0;
   *mxstrip = -999;
   for (;ptr<pend;) {
      if ( ptr->sn<SeedT || ptr->excl ) { /* does not pass trigger thresh.  */
	 ptr++;
	 continue;
      }
      /* Here we found a seed of a possible cluster */
      nright=1;
      nleft =0;
      nstrip=0;
      sn = ptr->sn;
      sg = ptr->signal;
      ph = ptr;              /* Save poiinter to seed */
      pl = ptr-1;
      pr = (ptr++)+1;
      ph->clus=1;
      /* Starts searching forward for strips passing inclusion threshold */
      while ( pr<pend && pr->sn>IncT ) {
	 nright++;
	 ptr++;
	 sn += pr->sn;
	 sg += pr->signal;
	 pr->clus=1;
	 pr++;
      }
      /* Search backward */
      while ( pl>=data && pl->sn>IncT ) {
	 nleft++;
	 sn += pl->sn;
	 sg += pr->signal;
	 pl->clus=1;
	 pl--;
      }
      nstrip = nright+nleft;
      *mxstrip = ( nclus > *mxstrip ? nstrip : *mxstrip );
      if ( sn > clsT ) {
	 nclus++;
	 pclst->ns  = nstrip;
	 pclst->fst = pl+1;
	 pclst->hig = ph;
	 pclst->sn  = sn;
	 pclst->sig = sg;	
	 /*  eta   */
	 if ( nstrip==1 )
	    pclst->eta=-0.1;
	 else {
	    if ( ph==pclst->fst ) {
	       pclst->eta = ph->signal/(ph->signal+(ph+1)->signal);
	    }
	    else if ( ph == (pclst->fst+pclst->ns) ) {
	       pclst->eta = ph->signal/(ph->signal+(ph-1)->signal);
	    }
	    else {
	       if ( (ph-1)->signal>(ph+1)->signal )
		  pclst->eta = ph->signal/(ph->signal+(ph-1)->signal);
	       else
		  pclst->eta = ph->signal/(ph->signal+(ph+1)->signal);
	    }
	 }	
	 pclst++;
      }
   }
   return (nclus);
}

 void Analog::GetPedestals()
{
register int i,j;
    for (i=0;i<3;i++) {
        printf("...Starting pass %d on the datan",i);
        pedes(i);
     }
    printf("...Got pedestals. Start the analysisn");
}


 void Analog::Reset_histograms(int how)
{
  register int i,j;
  switch (how)
    {
    case 1:
      for (j=0;j<Nreg;j++) 
	for(i=0;i<10;i++)
	 if(his[20*(j+1)+i]!=NULL)  his[20*(j+1)+i]->Reset();
        for(i=0;i<3;i++) prof[i]->Reset();
	his2[0]->Reset();
	break;
    default: 
   for(i=0;i<200;i++) 
     if(his[i]!=NULL) his[i]->Reset();
   for(i=0;i<5;i++) 
     if(his2[i]!=NULL) his2[i]->Reset();
   for(i=0;i<10;i++) 
     if(prof[i]!=NULL) prof[i]->Reset();
      break;
    }

}

 void Analog::DrawEvent(int ix)
{
  register int i;
  char name[100];

  ana(ix,1,0);
  
  for(i=0;i<4;i++)  
    {
      sprintf(name,"event view %d %d",ix,i);  
  if(chis[i]==NULL)
    chis[i]=new TH1F(name,name,NNchan,0.5,NNchan+0.5);
  else
    chis[i]->SetTitle(name);
    }


  double max=-1e4;
  double min=1e4;
	   
  for(i=0;i<NNchan;i++)
    { 
      if(!chan[i].excl)
	{
	  max=chan[i].signal>max?chan[i].signal:max;
	  min=chan[i].signal<min?chan[i].signal:min;
	}
      chis[0]->SetBinContent(i+1,chan[i].signal);
      chis[1]->SetBinContent(i+1,chan[i].sd);
      if(chan[i].clus) chis[2]->SetBinContent(i+1,chan[i].signal);
    }
    min=min>0?min*0.9:min*1.1;
  
    chis[0]->SetMaximum(max*1.1);
    chis[0]->SetMinimum(min);
    chis[0]->Draw(); 
    chis[1]->SetLineColor(2);
    chis[1]->Draw("SAME");      
    //    chis[2]->SetLineColor(4);
    //    chis[2]->Draw("SAME");      
//   for(i=1;i<4;i++) chis[i]->Draw("SAME");
  
}

void CLUster::Print()
{
  printf("No of strips in cluster: %d (%d,%d) man=%d n",ns,fst->id,fst->id+ns,hig->id);
  printf("Signal/Noise=%f, signal=%f, eta=%fn",sn,sig,eta);
  
}

void CHannel::Print()
{
  printf("id %d , excluded %d , cluster %dn",id,excl,clus);
  printf("raw %f , pedestal %f, noise %fn",raw,pd,sd);
  printf("sn %f , signal %f n",sn,signal);
  
}


ROOT page - Class index - Top of the page

This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.