//
//  evio2root - analysis program for data collected during the 3/2011
//		beam test of the Gluex tagger microscope.
//
// - based on micro_prot.cc, as, 21-mar-2011
// - based on fileAnalyzer.cc, ejw, 20-jun-2008
//
//  author: richard.t.jones at uconn.edu
//  date: April 2, 2011
//

#include <evioUtil.hxx>
#include <unistd.h>
#include <vector>
#include <math.h>

// root stuff
#include <TFile.h>
#include <TTree.h>
#include <TH1.h>
#include <TProfile.h>

#include "evio2root.h"

int evtCount = 0;

// prototypes
void analyzeEvent(evio::evioDOMTree &eventTree);
void analyzeBank(evio::evioDOMNodeP bankPtr);
void initRootTree(const char* filename);

Double_t sqr(Double_t x) { return x*x; }

// Define bit masks
// F1 bit masks
#define F1_RAW_DATA          0xFF000000

// Elliptical cuts in the correlation plots used to select
// events with primary pulse in each individual channel.
Cut2d *cut0 = new Cut2d(0,1,914,419,280,30,2.8,0);
Cut2d *cut1 = new Cut2d(1,3,794,503,257,30,20.3,0);
Cut2d *cut2 = new Cut2d(1,2,1018,1236,342,71,50,0);
Cut2d *cut3 = new Cut2d(3,4,1246,453,283,47,8.6,0);
Cut2d *cut4 = new Cut2d(3,4,385,1164,352,66,86,0);
Cut2d *cut12 = new Cut2d(12,13,957,490,229,18,0,0);
Cut2d *cut13 = new Cut2d(12,13,390,824,246,25,90,0);
Cut2d *cut14 = new Cut2d(12,14,390,874,246,25,90,0);

int main(int argc, char **argv)
{
  if (argc < 2) {
      std::cout << "Usage: evio2root <input_file.evio>" << std::endl;
      exit(1);
  }

  char filename[128];
  strcpy(filename, argv[1]);

  try {

    // open the output root file, using input filename with .root extension
    std::cout << "open output root file " << filename << std::endl;

    char fname[128],fname1[128];
    snprintf(fname1,128,"%s",filename);
    char *basename = fname1;
    if (strrchr(basename,'/')) {
       basename = strrchr(basename,'/');
       ++basename;
    }
    strtok(basename,".");
    sprintf(fname,"%s.root",basename);
    
    initRootTree(fname);
    
    // create evio file channel object using first arg as file name
    evio::evioFileChannel chan(argv[1]);
    
    // open the file
    std::cout << "open input evio data file " << argv[1] << std::endl;
    chan.open();
    
    // loop over events
    while (chan.read() && evtCount < 999999999) {
      evtCount++;

      if (evtCount%100==0) {
        std::cout << "fileAnalyzer read " << evtCount << " events" 
                  << std::endl;
      }
      
      // create event tree and analyze
      evio::evioDOMTree eventTree(chan);
      analyzeEvent(eventTree);
    }
    
    // reached eof, close file
    chan.close();
    
    ROOTfile->cd();
    ROOTfile->Write();
    ROOTfile->Close();
    delete ROOTfile;
    std::cout << std::endl << "\nClose ROOT file" << std::endl;

  } catch (evio::evioException *e) {
    std::cerr << std::endl << e->toString() << std::endl << std::endl;
    exit(EXIT_FAILURE);
  }
  
  // done
  std::cout << std::endl 
            << " ***evio2root done after " << evtCount << " events***"
            << std::endl << std::endl;
  exit(EXIT_SUCCESS);
}


void analyzeEvent(evio::evioDOMTree &eventTree)
{
  if (verbose) {
    std::cout << std::endl
              << "analyzing event " << evtCount 
              << std::endl << std::endl;
  }

  // reset statistics for this event
  for (int i=0; i<16; ++i) {
    amax[i] = 0;
    amin[i] = 9999;
    tmax[i] = 0;
    tmin[i] = 0;
  }
  tcount = 0;

  // get list of all banks in event
  evio::evioDOMNodeListP bankList = eventTree.getNodeList();

  // analyze each bank
  std::for_each(bankList->begin(),bankList->end(),analyzeBank);

  if (verbose) {
    std::cout << " Filling root File " << std::endl;
  }

  // make cuts based on pulse height
  if (cut0->Pass(amax)) {
    atree->Fill();
  }
  else if (cut1->Pass(amax)) {
    atree->Fill();
  }
  else if (cut2->Pass(amax)) {
    atree->Fill();
  }
  else if (cut3->Pass(amax)) {
    atree->Fill();
  }
  else if (cut4->Pass(amax)) {
    atree->Fill();
  }
  else if (cut12->Pass(amax)) {
    atree->Fill();
  }
  else if (cut13->Pass(amax)) {
    atree->Fill();
  }
  else if (cut14->Pass(amax)) {
    atree->Fill();
  }
}

void analyzeBank(evio::evioDOMNodeP bankPtr)
{

  // dump bank info for all banks
  //  std::cout << std::left << std::hex 
  //            << "bank content type:  0x"
  //            << std::setw(6) << bankPtr->getContentType() 
  //            << std::dec << "   tag:  " 
  //            << std::setw(6) << bankPtr->tag 
  //            << std::hex << "   num:  0x" 
  //            << std::setw(6) << (int)bankPtr->num 
  //            << std::dec << std::endl;

  // decode banks
  switch (bankPtr->tag) {
    
  case 0: // header bank
    //    std::cout << "found event bank" << std::endl;
    break;

  case 3:
    {
      const std::vector<int32_t> *vec = bankPtr->getVector<int32_t>();
      if (vec==NULL) {
        std::cerr << "?unable to get header bank vector" << std::endl;
        return;
      }
      //  std::cout << std::dec << "found header bank:  run number "
      //            << (*vec)[0] << "   Event number " << (*vec)[1]
      //            << std::endl;
      break;
    }
    
  case 17: // prestart event
    {
      const std::vector<uint32_t> *vec = bankPtr->getVector<uint32_t>();
      if (vec==NULL) {
        std::cerr << "?unable to get prestart event bank vector" << std::endl;
        return;
      }
      runNo = (*vec)[1];
      //  std::cout << std::dec << "found prestart bank:  run number "
      //            << (*vec)[1] << "   time " << (*vec)[0]
      //            << std::endl;
      break;
    }

  case 0xC000: // event id
    {
      const std::vector<uint32_t> *vec = bankPtr->getVector<uint32_t>();
      if (vec==NULL) {
        std::cerr << "?unable to get event id vector" << std::endl;
        return;
      }
        eventNo = (*vec)[0];
      //  std::cout << std::dec << "found event id bank:  event number "
      //            << (*vec)[0] << "   class " << (*vec)[1]
      //            << std::endl;
      break;
    }

  case 4: // fadc data bank
    {
      const std::vector<uint32_t> *vec = bankPtr->getVector<uint32_t>();
      if (vec==NULL) {
        std::cerr << "?unable to get fadc bank vector" << std::endl;
        return;
      }
      
      // unsigned int data_type =       
      // std::cout << "found fadc bank:  bank length " << vec->size() 
      //           << "    block length " << int( (*vec)[vec->size()-1] & 0x3fffff) 
      //           << "    block number " << int((*vec)[0]&0x3ff)
      //           << "    slot number "  << int(((*vec)[0]>>22)&0x1f)
      //           << std::endl;
      
      // Loop over data
      
      unsigned int data_new_type = 0;
      unsigned int data_type     = 0;
      int channel = -1;
      unsigned int time_stamp = 0;

      if (verbose) {
        std::cout << " Vector Size = " << vec->size() << std::endl;
      }

      for (unsigned i = 0; i < vec->size(); i++) {
        if ( (*vec)[i] & 0x80000000 ) {        /* data type defining word */
            data_type = ((*vec)[i] & 0x78000000) >> 27;
	    // data_type == 4  - window raw data
            data_new_type = 1;
        }
        else {
            data_new_type = 0;
        };
	
	//	std::cout << "SASCHA Data type = " << data_type << std::endl; 
	
        if (data_type == 4) {   // WINDOW RAW DATA
          if (data_new_type == 1) {
            unsigned int fadc_width = ((*vec)[i] & 0xFFF);
            channel    = ((*vec)[i] & 0x7800000) >> 23;
	    //  std::cout << " Data Type = " << data_type 
            //            << " Fadc_width  = " << fadc_width 
            //            << " Channel = " << channel << std::endl;
	    time_stamp = 1;
          }
          else {
            unsigned int valid_1 = 1;
            unsigned int valid_2 = 1;
            unsigned int adc_1 = ((*vec)[i] & 0x1FFF0000) >> 16;
	    
            if ((*vec)[i] & 0x20000000 ) {
              valid_1 = 0;
	    }
            unsigned int adc_2 = ((*vec)[i] & 0x1FFF);
            if ( (*vec)[i] & 0x2000 ) {
              valid_2 = 0;
	    }
	    
	    if ((channel == 0) && (time_stamp == 14)) {
	      //      std::cout << "CHANNEL = " << channel
              //                <<  " Time  stamp " << time_stamp
              //                <<  " Amplitude = " << adc_1 
              //                << "    " << adc_2 << std::endl;
	      ach1->Fill(float(adc_1));
	      ach2->Fill(float(adc_2));
	    }

            if (channel < 16) {
              if (adc_1 > amax[channel]) {
                amax[channel] = adc_1;
                tmax[channel] = time_stamp*2-2;
              }
              if (adc_1 < amin[channel]) {
                amin[channel] = adc_1;
                tmin[channel] = time_stamp*2-2;
              }
              if (adc_2 > amax[channel]) {
                amax[channel] = adc_2;
                tmax[channel] = time_stamp*2-1;
              }
              if (adc_2 < amin[channel]) {
                amin[channel] = adc_2;
                tmin[channel] = time_stamp*2-1;
              }
              if (evtCount <= 50) {
                ch[channel][evtCount-1]->Fill(float(time_stamp*2-2),float(adc_1),1.);
                ch[channel][evtCount-1]->Fill(float(time_stamp*2-1),float(adc_2),1.);
              }
            }

	    switch(channel){
	    case 0: 
	      ch0->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch0->Fill(float(time_stamp*2-1),float(adc_2),1.);
              trace0[time_stamp*2-2] = adc_1;
              trace0[time_stamp*2-1] = adc_2;
              break;
	    case 1:
	      ch1->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch1->Fill(float(time_stamp*2-1),float(adc_2),1.);
              trace1[time_stamp*2-2] = adc_1;
              trace1[time_stamp*2-1] = adc_2;
              break;
	    case 2:
	      ch2->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch2->Fill(float(time_stamp*2-1),float(adc_2),1.);
              trace2[time_stamp*2-2] = adc_1;
              trace2[time_stamp*2-1] = adc_2;
              break;
	    case 3:
	      ch3->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch3->Fill(float(time_stamp*2-1),float(adc_2),1.);
              trace3[time_stamp*2-2] = adc_1;
              trace3[time_stamp*2-1] = adc_2;
              break;
	    case 4:
	      ch4->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch4->Fill(float(time_stamp*2-1),float(adc_2),1.);
              trace4[time_stamp*2-2] = adc_1;
              trace4[time_stamp*2-1] = adc_2;
              break;
	    case 5:
	      ch5->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch5->Fill(float(time_stamp*2-1),float(adc_2),1.);
              break;
	    case 6:
	      ch6->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch6->Fill(float(time_stamp*2-1),float(adc_2),1.);
              break;
	    case 7:
	      ch7->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch7->Fill(float(time_stamp*2-1),float(adc_2),1.);
              break;
	    case 8:
	      ch8->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch8->Fill(float(time_stamp*2-1),float(adc_2),1.);
              break;
	    case 9:
	      ch9->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch9->Fill(float(time_stamp*2-1),float(adc_2),1.);
              break;
	    case 10:
	      ch10->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch10->Fill(float(time_stamp*2-1),float(adc_2),1.);
              break;
	    case 11:
	      ch11->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch11->Fill(float(time_stamp*2-1),float(adc_2),1.);
              break;
	    case 12:
	      ch12->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch12->Fill(float(time_stamp*2-1),float(adc_2),1.);
              trace12[time_stamp*2-2] = adc_1;
              trace12[time_stamp*2-1] = adc_2;
              break;
	    case 13:
	      ch13->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch13->Fill(float(time_stamp*2-1),float(adc_2),1.);
              trace13[time_stamp*2-2] = adc_1;
              trace13[time_stamp*2-1] = adc_2;
              break;
	    case 14:
	      ch14->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch14->Fill(float(time_stamp*2-1),float(adc_2),1.);
              trace14[time_stamp*2-2] = adc_1;
              trace14[time_stamp*2-1] = adc_2;
              break;
	    case 15:
	      ch15->Fill(float(time_stamp*2-2),float(adc_1),1.);
              ch15->Fill(float(time_stamp*2-1),float(adc_2),1.);
              trace15[time_stamp*2-2] = adc_1;
              trace15[time_stamp*2-1] = adc_2;
              break;
	    default:
	      break;
	    }

	    time_stamp++;

	    //	    if (channel == 10) {
	    //	      printf("%8X - RAW SAMPLES - valid = %d  adc = %d   valid = %d  adc = %d   channel = %d\n", 
	    //		     (*vec)[i], valid_1, adc_1, valid_2, adc_2, channel);	      	    	      
	    //	    }
	    //	    std::cout << " Time Stamp = " << time_stamp << std::endl;
            //      std::cout << std::hex << (*vec)[i] << std::endl;

          }
        }  // Data type 4, raw window data
	
	// Pulse Integral data
	if (data_type == 7) {
	  channel    = ((*vec)[i] & 0x7800000) >> 23;
	  unsigned int pulse_numb = ((*vec)[i] & 0x200000) >> 21;
	  unsigned int pulse_int  = ((*vec)[i] & 0x3FFFF);
	  
	  //  std::cout << " Channel = " << channel 
          //            << " Pulse number = " << pulse_numb
          //            << " Pulse int = " << pulse_int << std::endl;
	  // pulse_int = pulse_int - 13.4*12.;

	  switch (channel) {
	  case 0:
	    pc0->Fill(float(pulse_int));
            break;
	  case 1:
	    pc1->Fill(float(pulse_int));
            break;   
	  case 2:
	    pc2->Fill(float(pulse_int));
            break;
	  case 3:
	    pc3->Fill(float(pulse_int));
            break;
	  case 4:
	    pc4->Fill(float(pulse_int));
            break;
	  case 5:
	    pc5->Fill(float(pulse_int));
            break;
	  case 6:
	    pc6->Fill(float(pulse_int));
            break;
	  case 7:
	    pc7->Fill(float(pulse_int));
            break;
	  case 8:
	    pc8->Fill(float(pulse_int));
            break;
	  case 9:
	    pc9->Fill(float(pulse_int));
            break;
	  case 10:
	    pc10->Fill(float(pulse_int));
            break;
	  case 11:
	    pc11->Fill(float(pulse_int));
            break;
	  case 12:
	    pc12->Fill(float(pulse_int));
            break;
	  case 13:
	    pc13->Fill(float(pulse_int));
            break;
	  case 14:
	    pc14->Fill(float(pulse_int));
            break;
	  case 15:
	    pc15->Fill(float(pulse_int));
            break;
	  default:
	    break;	    
	  }
	}
      }
      
      break;
    }

  case 5: // F1TDC data
    {
      const std::vector<uint32_t> *vec = bankPtr->getVector<uint32_t>();
      if (vec==NULL) {
        std::cerr << "?unable to get tdc bank vector" << std::endl;
        return;
      }
      
      // Loop over data
      unsigned  int  data_new_type = 0;
      unsigned  int  data_type     = 0;
      int channel = -1;
      unsigned int time_stamp = 0;
      unsigned int F1_SLOT = 6;

      if (verbose) {
        std::cout << " -------------------" << std::endl;
        std::cout << " F1TDC Vector Size = " << vec->size() << std::endl;
        std::cout << " -------------------" << std::endl;
      }

      int time1 = -10;
      int time2 = -10;

      for (unsigned i = 0; i < vec->size(); i++) {
        if (verbose) {
	  std::cout  <<  "------------"   << std::endl;
	  std::cout  <<  std::hex << (*vec)[i] << std::dec << std::endl;
	  std::cout  <<  "------------"   << std::endl;
        }

	unsigned int slot_id = ((*vec)[i] & 0xF8000000) >> 27;
	unsigned int chip_res_locked = ((*vec)[i] & 0x4000000) >> 26;
	unsigned int overflow = ((*vec)[i] & 0x3000000) >> 24;
	unsigned int data_type = ((*vec)[i] & 0x800000) >> 23;

        if (verbose) {
	  std::cout << "   Slot id = " << slot_id 
                    << " chip_res_locked = " << chip_res_locked 
                    << "   Overflow = " << overflow 
                    << "  Data_type = " << data_type
                    << std::endl;
        }

	unsigned int f1_chip = 100;
	unsigned int f1_chan = 100;
	unsigned int f1_time = 0; 

	if (data_type == 0) {	  
	  unsigned int event_number = ((*vec)[i] & 0x3f0000) >> 16;
	  unsigned int trigger_time = ((*vec)[i] & 0xff80) >> 7;
	  f1_chip = ((*vec)[i] & 0x38) >> 3;
	  f1_chan = ((*vec)[i] & 0x7);  

          if (verbose) {
	    std::cout << std::endl;
	    std::cout << " Event number = " << event_number
                      << "  Trigger time = " << trigger_time
                      << " Chip = " << f1_chip
                      << "    Channel = " << f1_chan
                      << std::endl;
          }

	}
        else if (data_type == 1) { // Data
	   f1_chip = ((*vec)[i] & 0x380000) >> 19;
	   f1_chan = ((*vec)[i] & 0x70000) >> 16;
	   f1_time = ((*vec)[i] & 0xffff);
           if (verbose) {
	     std::cout << " Chip = " << f1_chip
                       << "      Channel = " << f1_chan
                       << "      Time = " << f1_time
                       << std::endl;
           }
	   if ((f1_chip == 0) && ((f1_chan == 2) || (f1_chan == 3))) {
             time1 = f1_time;
           }
	   if ((f1_chip == 3) && ((f1_chan == 6) || (f1_chan == 7))) {
             time2 = f1_time;
           }
           tchip[tcount] = f1_chip;
           tchan[tcount] = f1_chan;
           ttime[tcount] = f1_time;
           tcount++;
	}

	if (slot_id != F1_SLOT) {
          // std::cerr << " Wrong slot number " << slot_id << std::endl;
        }
      }

      if (verbose) {
        std::cout << " Time difference = " << time2 - time1 << std::endl;
      }
      time_diff->Fill(float(time2-time1));

    }

  default:
      break;
  }
}

void initRootTree(const char* fname)
{
  ROOTfile = new TFile(fname,"RECREATE","new");
  std::cout << "Opened ROOT file " << fname << " ..." << std::endl;
  std::cout << "Create Root tree ..." << std::endl;
  ROOTfile->cd();

  atree = new TTree("FadcStuff","Fadc stuff");
  atree->Branch("runNo", &runNo, "runNo/I");
  atree->Branch("eventNo", &eventNo, "eventNo/I");
  atree->Branch("amax", amax, "amax[16]/S");
  atree->Branch("amin", amin, "amin[16]/S");
  atree->Branch("tmax", tmax, "tmax[16]/S");
  atree->Branch("tmin", tmin, "tmin[16]/S");
  atree->Branch("trace0", trace0, "trace0[125]/S");
  atree->Branch("trace1", trace1, "trace1[125]/S");
  atree->Branch("trace2", trace2, "trace2[125]/S");
  atree->Branch("trace3", trace3, "trace3[125]/S");
  atree->Branch("trace4", trace4, "trace4[125]/S");
  atree->Branch("trace12", trace12, "trace12[125]/S");
  atree->Branch("trace13", trace13, "trace13[125]/S");
  atree->Branch("trace14", trace14, "trace14[125]/S");
  atree->Branch("trace15", trace15, "trace15[125]/S");
  atree->Branch("tcount", &tcount, "tcount/I");
  atree->Branch("tchip", tchip, "tchip[30]/S");
  atree->Branch("tchan", tchan, "tchan[30]/S");
  atree->Branch("ttime", ttime, "ttime[30]/S");

  int window_size = 125;

  ch0 = new TProfile("ch 0","ch 0",window_size,-0.5,window_size-0.5,-10.,4096.);
  ch1 = new TProfile("ch 1","ch 1",window_size,-0.5,window_size-0.5,-10.,4096.);
  ch2 = new TProfile("ch 2","ch 2",window_size,-0.5,window_size-0.5,-10.,4096.);
  ch3 = new TProfile("ch 3","ch 3",window_size,-0.5,window_size-0.5,-10.,4096.);
  ch4 = new TProfile("ch 4","ch 4",window_size,-0.5,window_size-0.5,-10.,4096.);
  ch5 = new TProfile("ch 5","ch 5",window_size,-0.5,window_size-0.5,-10.,4096.);
  ch6 = new TProfile("ch 6","ch 6",window_size,-0.5,window_size-0.5,-10.,4096.);
  ch7 = new TProfile("ch 7","ch 7",window_size,-0.5,window_size-0.5,-10.,4096.);
  ch8 = new TProfile("ch 8","ch 8",window_size,-0.5,window_size-0.5,-10.,4096.);
  ch9 = new TProfile("ch 9","ch 9",window_size,-0.5,window_size-0.5,-10.,4096.);
  ch10 = new TProfile("ch 10","ch 10",window_size,-0.5,window_size-0.5,-10.,4096.);
  ch11 = new TProfile("ch 11","ch 11",window_size,-0.5,window_size-0.5,-10.,4096.);
  ch12 = new TProfile("ch 12","ch 12",window_size,-0.5,window_size-0.5,-10.,4096.);
  ch13 = new TProfile("ch 13","ch 13",window_size,-0.5,window_size-0.5,-10.,4096.);
  ch14 = new TProfile("ch 14","ch 14",window_size,-0.5,window_size-0.5,-10.,4096.);
  ch15 = new TProfile("ch 15","ch 15",window_size,-0.5,window_size-0.5,-10.,4096.);


  pc0 = new TH1F("pch 0","pch 0",3400,-0.5,3399.5);
  pc1 = new TH1F("pch 1","pch 1",3400,-0.5,3399.5);
  pc2 = new TH1F("pch 2","pch 2",3400,-0.5,3399.5);
  pc3 = new TH1F("pch 3","pch 3",3400,-0.5,3399.5);
  pc4 = new TH1F("pch 4","pch 4",3400,-0.5,3399.5);
  pc5 = new TH1F("pch 5","pch 5",3400,-0.5,3399.5);
  pc6 = new TH1F("pch 6","pch 6",3400,-0.5,3399.5);
  pc7 = new TH1F("pch 7","pch 7",3400,-0.5,3399.5);
  pc8 = new TH1F("pch 8","pch 8",3400,-0.5,3399.5);
  pc9 = new TH1F("pch 9","pch 9",3400,-0.5,3399.5);
  pc10 = new TH1F("pch 10","pch 10",3400,-0.5,3399.5);
  pc11 = new TH1F("pch 11","pch 11",3400,-0.5,3399.5);
  pc12 = new TH1F("pch 12","pch 12",3400,-0.5,3399.5);
  pc13 = new TH1F("pch 13","pch 13",3400,-0.5,3399.5);
  //  pc14   =  new TH1F("pch 14","pch 14",10000,-0.5,9999.5);
  pc14 = new TH1F("pch 14","pch 14",40000,14999.5,54999.5);
  pc15 = new TH1F("pch 15","pch 15",300,-0.5,299.5);

  ach1 = new TH1F("ach 1","ach 1",1024,-0.5,1023.5);
  ach2 = new TH1F("ach 2","ach 2",1024,-0.5,1023.5);

  time_diff = new TH1F("Dt","Dt",100,120.,220.);

  for (int event=0; event<50; ++event) {
    for (int chan=0; chan<16; ++chan) {
      char name[10];
      sprintf(name,"c%de%d",chan,event);
      ch[chan][event] = new TProfile(name,name,window_size,-0.5,window_size-0.5,-10.,4096.);
    }
  }
}

Cut2d::Cut2d(Int_t c0, Int_t c1,
             Double_t x0, Double_t y0,
             Double_t amajor, Double_t aminor,
             Double_t theta, Int_t veto) 
{
   fChan[0] = c0;
   fChan[1] = c1;
   fOrigin[0] = x0;
   fOrigin[1] = y0;
   fAxis[0] = amajor;
   fAxis[1] = aminor;
   fThetaR = theta*3.1415926/180;
   fReject = veto;
}

Int_t Cut2d::Pass(short int *alist) const
{
   Double_t x = alist[fChan[0]]-fOrigin[0];
   Double_t y = alist[fChan[1]]-fOrigin[1];
   Double_t cost = cos(fThetaR);
   Double_t sint = sin(fThetaR);
   Double_t u = x*cost+y*sint;
   Double_t v = y*cost-x*sint;
   Double_t dist = sqrt(sqr(u/fAxis[0])+sqr(v/fAxis[1]));
   return (fReject)? (dist > 1) : (dist <= 1);
}