void calcthr()
{
   // This version modified 26-Nov-2003 by JPSullivan
   // It now does two passes in calculating the mean and rms for each channel,
   // the first pass uses everything. The second pass only keeps events with 
   // adc values within "close*sigma" were close is a variable currently set
   // to 3.0 and sigma is the rms width from the first pass. This is intended
   // to keep events near the preamp reset out of the calculation of the pedestal
   // location.
   //
   // This routine does all of its calculations based on an input file called
   // check.root, which is created (before calling this routine) by calling
   // check.C

   printf ( " calcthr.C enter: setup threshold files\n");
   gROOT->Reset();
   TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("check.root");
   if (!f) 
   {
      f = new TFile("check.root");
   }
   TTree *ntp = (TTree*)gDirectory->Get("ntp");
   
   //Declaration of leaves types
   Float_t         evtseq;
   Float_t         packetid;
   Float_t         user7;
   Float_t         amu0;
   Float_t         amu1;
   Float_t         channel;
   Float_t         adc;
   
   int ipid, ich;
   const int packetstart = 2001;
   const int packetstop = 2120;
   const int npackets = packetstop - packetstart + 1;
   const int NCHANNELS = 256;
   const float notvalid = -99;
   
   float thrmean[npackets][NCHANNELS];
   float thrsigma[npackets][NCHANNELS];
   int threntries[npackets][NCHANNELS];
   float thrsum[npackets][NCHANNELS];
   float thrsum2[npackets][NCHANNELS];

   for (int ipacketid = packetstart; ipacketid<=packetstop; ipacketid++) 
   {
      ipid = ipacketid - packetstart;
      for (ich = 0; ich<NCHANNELS; ich++)
      {
         thrmean[ipid][ich]    = 0.0;
         thrsigma[ipid][ich]   = 0.0;
         threntries[ipid][ich] = 0;
         thrsum[ipid][ich]     = 0.0;
         thrsum2[ipid][ich]    = 0.0;
      }
   }
   
   //Set branch addresses
   ntp->SetBranchAddress("evtseq",&evtseq);
   ntp->SetBranchAddress("packetid",&packetid);
   ntp->SetBranchAddress("user7",&user7);
   ntp->SetBranchAddress("amu0",&amu0);
   ntp->SetBranchAddress("amu1",&amu1);
   ntp->SetBranchAddress("channel",&channel);
   ntp->SetBranchAddress("adc",&adc);

   Int_t nentries = ntp->GetEntries();
   
   // add up the sum of the adc and the adc squared
   Int_t nbytes = 0;
   float pid;
   for (Int_t i=0; i<nentries;i++) 
   {
      nbytes += ntp->GetEntry(i);
      pid = packetid - packetstart;
      if (pid >= 0 && pid<npackets)
      {
         ipid = (int) pid;
         ich = (int) channel;

         thrsum[ipid][ich]  += adc;
         thrsum2[ipid][ich] += adc*adc;
         threntries[ipid][ich]++;
      }
   }


   float rms;
   // calculate the mean and sigma on a first pass -- without any cuts
   for (int ipacketid = packetstart; ipacketid<=packetstop; ipacketid++) 
   {
      ipid = ipacketid - packetstart;
      for (ich = 0; ich<NCHANNELS; ich++)
      {
         if ( threntries[ipid][ich] > 0 )
         {
            thrmean[ipid][ich] = thrsum[ipid][ich]/(float)threntries[ipid][ich];
            rms = thrsum2[ipid][ich]/(float)threntries[ipid][ich]
		    - thrmean[ipid][ich]*thrmean[ipid][ich];
            if ( rms > 0. )
            {
               thrsigma[ipid][ich] = sqrt(rms);
            }
            else
            {
               thrsigma[ipid][ich] = notvalid;
            }
         }
         else
         {
            thrmean[ipid][ich] = notvalid;
            thrsigma[ipid][ich] = notvalid;
         }
      }
   }

   printf ( " calcthr.C: end first pass\n");

   // reset the numbers used in the calculation of the mean and sigma
   for (int ipacketid = packetstart; ipacketid<=packetstop; ipacketid++) 
   {
      ipid = ipacketid - packetstart;
      for (ich = 0; ich<NCHANNELS; ich++)
      {
         threntries[ipid][ich] = 0;
         thrsum[ipid][ich]=0.0;
         thrsum2[ipid][ich]=0.0;
      }
   }
   // Now pass through the data again, but only keep the values which
   // are close to the mean value (mainly used to reject skewed means caused
   // by events near the preamp reset)
   float close = 3.;  // close*sigma defines "close to" the mean
   nbytes = 0;
   for (Int_t i=0; i<nentries;i++) 
   {
      nbytes += ntp->GetEntry(i);
      pid = packetid - packetstart;
      if (pid >= 0 && pid<npackets)
      {
         ipid = (int) pid;
         ich = (int) channel;
         if ( fabs(adc - thrmean[ipid][ich]) < close*thrsigma[ipid][ich] )
         {
            thrsum[ipid][ich] += adc;
            thrsum2[ipid][ich] += adc*adc;
            threntries[ipid][ich]++;
         }
      }
   }

   // calculate the mean and sigma -- with cuts
   float old_mean, old_rms;
   for (int ipacketid = packetstart; ipacketid<=packetstop; ipacketid++) 
   {
      ipid = ipacketid - packetstart;
      printf ( " ipid=%d\n",ipid);
      for (ich = 0; ich<NCHANNELS; ich++)
      {
         old_mean = thrmean[ipid][ich];
         old_rms   = thrsigma[ipid][ich];
         if ( threntries[ipid][ich] > 0 )
         {
            thrmean[ipid][ich] = thrsum[ipid][ich]/(float)threntries[ipid][ich];
            rms = thrsum2[ipid][ich]/(float)threntries[ipid][ich]
		    - thrmean[ipid][ich]*thrmean[ipid][ich];
            if ( rms > 0. )
            {
               thrsigma[ipid][ich] = sqrt(rms);
            }
            else
            {
               thrsigma[ipid][ich] = notvalid;
            }
         }
         else
         {
            thrmean[ipid][ich] = notvalid;
            thrsigma[ipid][ich] = notvalid;
         }
         if ( fabs(old_mean-thrmean[ipid][ich])>thrsigma[ipid][ich] && thrsigma[ipid][ich]>0. )
         {
             // print a warning if there is a big shift in the pedestal between
             // the first and 2nd passes, this implies something strange going on.
             printf ( " calcthr.C: more than 1 sigma change in 2nd pass ped calc" );
             printf ( " ipacketid=%d ich=%d 1st pass mean=%f rms=%f 2nd pass mean=%f rms=%f\n",
                      ipid+1, ich, old_mean, old_rms, thrmean[ipid][ich], thrsigma[ipid][ich]);
         }
      }
   }

   printf ( " calcthr.C: end second pass\n");

   // create the report
   char outname[20];
   float cut; 
   // const float sigcut = 3.0; // until March 7, 2003; and back again
   // const float sigcut = 2.0; // after March 7, 2003
   const float sigcut    = 2.0; // after Nov 25, 2003 (aka start run 4)
   const float maxabove  = 10.0; // the cut will not be more than this number of channels
                                 // above the pedestal in any case 
   // use this for AuAu: const int chansave    = 8;   // save every 8 channels (not 0-supp)
   const int chansave    = 16;   // save every 16 channels (not 0-supp)

   FILE* fsum=fopen("mvdsummary.txt","w");

   for (int ipacketid = packetstart; ipacketid<=packetstop; ipacketid++) 
   {
      ipid = ipacketid - packetstart;
      sprintf(outname,      "mvd%d.thresh",ipacketid);   
      FILE* fout=fopen(outname,"w");

      for (ich = 0; ich<NCHANNELS; ich++)
      {
         
         if ( sigcut*thrsigma[ipid][ich] < maxabove ) 
         { // "normal"
            cut = thrmean[ipid][ich] + sigcut*thrsigma[ipid][ich];  
         }
         else 
         {
            cut = thrmean[ipid][ich] + maxabove;   
            // cout  << " limit cut to 10 chan above pedestal" << endl;
         }
         
         // adjustments
         if ( thrsigma[ipid][ich] > 30 ) 
         { // too noisy - shouldn't these be masked out?
            cut = 3;
         }
         if ( cut < 0 ) 
         {
            cut = 1022;  // bad fits, or invalid data 
         }
         if ( cut > 1024) 
         { // too high
            cut = 4;
         }               
         // leave some channels not zero-suppressed
         if ( ich%chansave == 0 )  
         {
            cut = 1;
         }

         fprintf(fout,"%d\n",int(cut+0.5));
         fprintf(fsum,"%d %3d %6.2f %6.2f %d\n",ipacketid,ich,
                     thrmean[ipid][ich],thrsigma[ipid][ich], int(cut+0.5));
      }
      fclose(fout);
   }
   fclose(fsum);
   printf ( " calcthr.C: output files created\n");
   return;
}