IceRoot.cxx

Go to the documentation of this file.

 
#include "IceRoot.h"
#include "Riostream.h"
 
ClassImp(IceRoot); // Class implementation to enable ROOT I/O
 
IceRoot::IceRoot(const char* name,const char* title) : NcJob(name,title)
{
 
 fSplit=0;
 fBsize=32000;
 fMaxevt=-1;
 fPrintfreq=1;
 fInfiles=0;
 fOutfile=0;
 fCalfile=0;
 fJEBTDaq=0;
}

IceRoot::~IceRoot()
{
 
 if (fInfiles)
 {
  delete fInfiles;
  fInfiles=0;
 }
 
 if (fCalfile)
 {
  delete fCalfile;
  fCalfile=0;
 }
}

void IceRoot::SetMaxEvents(Int_t n)
{
 
 fMaxevt=n;
}

void IceRoot::SetPrintFreq(Int_t f)
{
 
 if (f>=0) fPrintfreq=f;
}

void IceRoot::SetSplitLevel(Int_t split)
{
 
 if (split>=0) fSplit=split;
}

void IceRoot::SetBufferSize(Int_t bsize)
{
 
 if (bsize>=0) fBsize=bsize;
}

void IceRoot::SetInputFile(TString name)
{
 
 if (!fInfiles)
 {
  fInfiles=new TObjArray();
  fInfiles->SetOwner();
 }
 fInfiles->Clear();
 
 TObjString* s=new TObjString();
 s->SetString(name);
 fInfiles->Add(s);
}

void IceRoot::AddInputFile(TString name)
{
 
 if (!fInfiles)
 {
  fInfiles=new TObjArray();
  fInfiles->SetOwner();
 }
 
 TObjString* s=new TObjString();
 s->SetString(name);
 fInfiles->Add(s);
}

void IceRoot::SetOutputFile(TFile* ofile)
{
 
 if (fOutfile) delete fOutfile;
 fOutfile=ofile;
}

void IceRoot::SetOutputFile(TString name)
{
 
 if (fOutfile) delete fOutfile;
 fOutfile=new TFile(name.Data(),"RECREATE","Simple Root data in IceEvent structure");
}

TFile* IceRoot::GetOutputFile()
{
 
 return fOutfile;
}

void IceRoot::SetCalibFile(TString name)
{
 
 if (fCalfile) delete fCalfile;
 fCalfile=new TFile(name.Data());
 
 if(fCalfile){
  fJEBTDaq=(NcObjMatrix*)fCalfile->Get("JEBTDaq-OMDBASE");
 }
 if(!fJEBTDaq){
  cout << "*IceRoot* Warning: no calibration available for TWR. TWR waveforms cannot be processed." << endl;
 }
 
}

void IceRoot::SetOMdbase(NcObjMatrix* omdb)
{
 
 fJEBTDaq=omdb;
 if(!fJEBTDaq){
  cout << "*IceRoot* Warning: no calibration available for TWR. TWR waveforms cannot be processed." << endl;
 }
 
}

void IceRoot::Exec(Option_t* opt)
{
 
 if (!fInfiles)
 {
  cout << " *IceRoot Exec* No data input file(s) specified." << endl;
  return;
 }
 
 Int_t ninfiles=fInfiles->GetEntries();
 if (!ninfiles)
 {
  cout << " *IceRoot Exec* No data input file(s) specified." << endl;
  return;
 }
 
 // Warning in case no calibration DB is specified
 if(!fJEBTDaq){
  cout << "*IceRoot* Warning: no calibration available for TWR. TWR waveforms cannot be processed." << endl;
 }
 
 // Create output tree if necessary
 TTree* otree=0;
 if (fOutfile)
 {
  otree=new TTree("T","Simple Root data converted to IceEvent structures");
  otree->SetDirectory(fOutfile);
 }
 
 // Create IceEvent structure
 IceEvent* evt=new IceEvent();
 evt->SetTrackCopy(1);
 evt->SetDevCopy(1);
 
 // Branch in the tree for the event structure
 if (otree) otree->Branch("IceEvent","IceEvent",&evt,fBsize,fSplit); 
 
 // Initialise the job working environment
 SetMainObject(evt);
 
 // Some output for the user's convenience
 TString inputfile;
 cout << " ***" << endl;
 cout << " *** Start processing of job " << GetName() << " ***" << endl;
 cout << " ***" << endl;
 for (Int_t i=0; i<ninfiles; i++)
 {
  TObjString* sx=(TObjString*)fInfiles->At(i);
  if (!sx) continue;
  inputfile=sx->GetString(); 
  cout << " Simple Root data input file : " << inputfile.Data() << endl;
 }
 cout << " Maximum number of events to be processed : " << fMaxevt << endl;
 cout << " Print frequency : " << fPrintfreq << endl;
 if (fOutfile)
 {
  cout << " ROOT output file : " << fOutfile->GetName() << endl;
  cout << " Output characteristics : splitlevel = " << fSplit << " buffersize = " << fBsize << endl;
 }
 
 ListEnvironment();
 
 // Set DAQ device info
 NcDevice daq;
 daq.SetNameTitle("Daq","Daq system used");
 
 // Trigger device
 NcDevice trig;
 trig.SetNameTitle("Trigger","Amanda/IceCube event triggers");
 NcSignal s;
 Char_t trigname[100];
 
 // MC info device
 NcDevice mcinfo;
 mcinfo.SetNameTitle("MCInfo","The MC production info");
 
 // Variables in the simple ROOT tree
 Int_t     pretrig=0;
 Int_t     trignr=0;
 Double_t* trigtime=0;
 Double_t* triglength=0;
 Char_t    trigsourceID[100][100];
 Char_t    trigtypeID[100][100];
 Char_t    trigsubtypeID[100][100];
 Double_t* mcweightitemvalue=0;
 Char_t    mcweightitemname[100][100];
 Int_t     eventtimemjd=0;
 Int_t     eventtimemjds=0;
 Int_t     eventtimemjdns=0;
 Int_t     eventid=0;
 Int_t     runid=0;
 Int_t     String=0;
 Int_t     OM=0;
 Short_t   waveformtype=0; // 1 = TWR, 2 = non-merged ATWD (not supported), 3 = InIce merged ATWD, 4 = InIce FADC, 5 = IceTop merged ATWD, 6 = IceTop FADC
 Short_t   twreventnr=0;
 UInt_t    baseline=0;
 Int_t     numberbinstwr=0;
 Int_t     startbin=0;
 UInt_t*   twrwaveform=0;
 Int_t     numberbinswaveform=0;
 Double_t  starttimewaveform=0;
 Double_t  binwidth=0;
 Short_t   source=0; // 0 = ATWD, 10 = FADC, 20 = TWR elec, 30 = TWR opt, 40 = etc
 Double_t* waveform=0;
 Int_t     ntracks=0;
 Double_t  paraboloidsigmaA=999;
 Double_t  paraboloidsigmaI=999;
 Double_t  rloglA=999;
 Double_t  rloglI=999;
 Double_t* trackx=0;
 Double_t* tracky=0;
 Double_t* trackz=0;
 Double_t* trackzenith=0;
 Double_t* trackazimuth=0;
 Double_t* trackenergy=0;
 Char_t    tracktype[1000][100];
 
 // Some other variables
 Int_t nevt=0;
 Int_t lastevent=0;
 Int_t omid=0;
 Double_t starttime=0;
 Int_t extstop=0;
 Float_t twrtime=10240; // Number of bins in TWR
 Int_t twrbuffer=0;
 TPRegexp suffix("_[0-9]+$"); // suffix appended to track name by Walnut
 
 TString hname;
 TH1F histo;
 IceAOM aom;
 IceIDOM idom;
 IceTDOM tdom;
 IceGOM* om=0;
 IceGOM* calom=0;
 NcTrack t;
 Double_t vec[3];
 NcPosition r;
 Nc3Vector p;
 NcSample sample;
 
 Float_t pi=acos(-1.);
 
 TFile* input=0;
 
 Double_t* twreventtimes=0;
 Int_t nrtwrevents=0;
 Int_t* twreventorder=0;
 
 // Get global TWR time offsets from DB
 Float_t globalt0=0;
 Float_t twri3offset=0;
 if(fJEBTDaq){
  for(Int_t row=1; row<=fJEBTDaq->GetMaxRow(); row++){
   calom=(IceGOM*)fJEBTDaq->GetObject(row,1);
   if(calom){
    globalt0=calom->GetSignal("GLOBALT0");
    twri3offset=calom->GetSignal("TWRI3OFFSET");
    break;
   }
  }
 }
 
 for (Int_t ifile=0; ifile<ninfiles; ifile++)
 {
  TObjString* sx=(TObjString*)fInfiles->At(ifile);
  if (!sx) continue;
 
  inputfile=sx->GetString(); 
  if (inputfile=="") continue;
 
  // Open the simple Root data input file
  input=TFile::Open(inputfile);
 
  if (!input)
  {
   cout << " *IceRoot Exec* No input file found with name : " << inputfile.Data() << endl;
   continue;
  }
 
  // Get simple Root tree
  fTree=(TTree*)input->Get("T");
 
  fTree->SetBranchAddress("pretrig",&pretrig);
  fTree->SetBranchAddress("trignr",&trignr);
  fTree->SetBranchAddress("eventtimemjd",&eventtimemjd);
  fTree->SetBranchAddress("eventtimemjds",&eventtimemjds);
  fTree->SetBranchAddress("eventtimemjdns",&eventtimemjdns);
  fTree->SetBranchAddress("eventid",&eventid);
  fTree->SetBranchAddress("runid",&runid);
  fTree->SetBranchAddress("String",&String);
  fTree->SetBranchAddress("OM",&OM);
  fTree->SetBranchAddress("waveformtype",&waveformtype);
  fTree->SetBranchAddress("twreventnr",&twreventnr);
  fTree->SetBranchAddress("baseline",&baseline);
  fTree->SetBranchAddress("numberbinstwr",&numberbinstwr);
  fTree->SetBranchAddress("startbin",&startbin);
  fTree->SetBranchAddress("numberbinswaveform",&numberbinswaveform);
  fTree->SetBranchAddress("starttimewaveform",&starttimewaveform);
  if(fTree->GetBranch("binwidth")) fTree->SetBranchAddress("binwidth",&binwidth);
  fTree->SetBranchAddress("source",&source);
  fTree->SetBranchAddress("ntracks",&ntracks);
 
  Int_t ntrig=(Int_t)fTree->GetLeaf("trignr")->GetMaximum();
  trigtime=new Double_t[ntrig+1];
  fTree->SetBranchAddress("trigtime",trigtime);
  triglength=new Double_t[ntrig+1];
  fTree->SetBranchAddress("triglength",triglength);
  fTree->SetBranchAddress("trigsourceID",trigsourceID);
  fTree->SetBranchAddress("trigtypeID",trigtypeID);
  fTree->SetBranchAddress("trigsubtypeID",trigsubtypeID);
 
  Int_t nmcinfo=(Int_t)fTree->GetLeaf("Nmcinfo")->GetMaximum();
  mcweightitemvalue=new Double_t[nmcinfo+1];
  fTree->SetBranchAddress("mcweightitemvalue",mcweightitemvalue);
  fTree->SetBranchAddress("mcweightitemname",mcweightitemname);
 
  Int_t nmaxbinstwr=(Int_t)fTree->GetLeaf("numberbinstwr")->GetMaximum();
  twrwaveform=new UInt_t[nmaxbinstwr+1];
  fTree->SetBranchAddress("twrwaveform",twrwaveform);
 
  Int_t nmaxbinswaveform=(Int_t)fTree->GetLeaf("numberbinswaveform")->GetMaximum();
  waveform=new Double_t[nmaxbinswaveform+1];
  fTree->SetBranchAddress("waveform",waveform);
 
  fTree->SetBranchAddress("paraboloidsigmaA",&paraboloidsigmaA);
  fTree->SetBranchAddress("paraboloidsigmaI",&paraboloidsigmaI);
  fTree->SetBranchAddress("rloglA",&rloglA);
  fTree->SetBranchAddress("rloglI",&rloglI);
 
  Int_t nmaxtracks=(Int_t)fTree->GetLeaf("ntracks")->GetMaximum();
  trackx=new Double_t[nmaxtracks];
  tracky=new Double_t[nmaxtracks];
  trackz=new Double_t[nmaxtracks];
  trackzenith=new Double_t[nmaxtracks];
  trackazimuth=new Double_t[nmaxtracks];
  trackenergy=new Double_t[nmaxtracks];
  fTree->SetBranchAddress("trackx",trackx);
  fTree->SetBranchAddress("tracky",tracky);
  fTree->SetBranchAddress("trackz",trackz);
  fTree->SetBranchAddress("trackzenith",trackzenith);
  fTree->SetBranchAddress("trackazimuth",trackazimuth);
  fTree->SetBranchAddress("trackenergy",trackenergy);
  fTree->SetBranchAddress("tracktype",tracktype);
 
  // Prepare for loop over entries
  lastevent=0;
 
  // Loop over waveforms in input tree
  for(Int_t ientry=0; ientry<fTree->GetEntries(); ientry++)
  {
   fTree->GetEntry(ientry);
   // If new event
   if(eventid!=lastevent)
   {
    // Write old event to tree (if it contains data)
    if(evt->GetDevices("IceGOM"))
    {
     // Invoke all available sub-tasks (if any) and write event to tree
     CleanTasks();
     ExecuteTasks(opt);
     if(otree) otree->Fill();
     if (fPrintfreq) { if (!(nevt%fPrintfreq)) evt->HeaderData(); }
     // Update event counter
     nevt++;
     if (fMaxevt>-1 && nevt>=fMaxevt) 
     {
      evt->Reset();
      break;
     }
    }
 
    // Start new event
    evt->Reset();
    evt->SetRunNumber(runid);
    evt->SetEventNumber(eventid);
    evt->SetMJD(eventtimemjd,eventtimemjds,eventtimemjdns,0);
 
    // Daq: JEB for 2007 and later, TWR for 2006 and earlier
    daq.Reset(1);
    if(evt->GetDate()/10000 >= 2007){
     daq.AddNamedSlot("JEB");
     daq.SetSignal(1,"JEB");
    }
    else {
     daq.AddNamedSlot("TWR");
     daq.SetSignal(1,"TWR");
    }
    evt->AddDevice(daq);
 
    // Loop over all the tracks and add them to the current event
    Int_t nrecotracks=0;
    Int_t nmctracks=0;
    for (Int_t itrack=0; itrack<ntracks; itrack++)
    {
     t.Reset();
 
     TString tracktypestring(tracktype[itrack]);
     // Monte Carlo track
     if(tracktypestring.Index("MC")>=0 || tracktypestring.Index("nu_mu")>=0){
      nmctracks++;
      t.SetId(-nmctracks);
     }
     // Reco track
     if(tracktypestring.Index("MC")<0 && tracktypestring.Index("nu_mu")<0){
      nrecotracks++;
      t.SetId(nrecotracks);
     }
     // Track ID and name
     Int_t suffixpos=tracktypestring.Index(suffix);
     if(suffixpos>0){
      tracktypestring.Remove(suffixpos,tracktypestring.Length()-suffixpos);
     }
     t.SetName(tracktypestring);
     t.SetTitle(tracktypestring);
 
     // Beginpoint of the track
     vec[0]=trackx[itrack];
     vec[1]=tracky[itrack];
     vec[2]=trackz[itrack];
     r.SetPosition(vec,"car");
     t.SetBeginPoint(r);
 
     // Momentum in GeV/c
     vec[0]=trackenergy[itrack];
     if(vec[0]==0 || vec[0]!=vec[0]) vec[0]=1; // Energy unknown: set default value of 1 GeV
     vec[1]=pi-trackzenith[itrack];
     vec[2]=trackazimuth[itrack]+pi;
     if(vec[2]>=2*pi) vec[2]-=2*pi;
     p.SetVector(vec,"sph");
     t.Set3Momentum(p);
 
     if(tracktypestring.Index("paraboloidA")>=0){
       t.SetCharge(paraboloidsigmaA);
     }else if(tracktypestring.Index("paraboloidI")>=0){
       t.SetCharge(paraboloidsigmaI);
     }
     if(tracktypestring.Index("gulliverA")>=0){
       t.SetCharge(rloglA);
     }else if(tracktypestring.Index("gulliverI")>=0){
       t.SetCharge(rloglI);
     }
 
     // Check for unreconstructed tracks (NaN values) and reset track if necessary
     if(trackx[itrack]!=trackx[itrack] || tracky[itrack]!=tracky[itrack] || trackz[itrack]!=trackz[itrack] ||
        trackzenith[itrack]!=trackzenith[itrack] || trackazimuth[itrack]!=trackazimuth[itrack]){
      t.Reset();
     } else {
      // Add track to event
      evt->AddTrack(t);
     }
    }
 
    // Store trigger information
    if (twreventtimes) delete[] twreventtimes;
    twreventtimes=new Double_t[trignr];
    if (twreventorder) delete[] twreventorder;
    twreventorder=new Int_t[trignr];
    nrtwrevents=0;
    trig.Reset(1);
    // Loop over triggers
    for(Int_t itr=0; itr<trignr; itr++){
     // For TWR events:
     if(TString(trigsourceID[itr])=="AMANDA_TWR_DAQ"){
      // Add global time offset if there are no MC tracks (i.e. if the data are real data)
      if(!nmctracks) trigtime[itr]-=globalt0+twri3offset;
      // Check if this trigger belongs to a new event
      Int_t newevent=1;
      for(Int_t itwrevent=0; itwrevent<nrtwrevents; itwrevent++){
       if(fabs(trigtime[itr]-twreventtimes[itwrevent])<1500){
        newevent=0;
        break;
       }
      }
      // Add new TWR event time
      if(newevent){
       twreventtimes[nrtwrevents]=trigtime[itr];
       nrtwrevents++;
      }
     }
 
     // Add this trigger to trigger device
     s.Reset(1);
     sprintf(trigname,"%s/%s",trigsourceID[itr],trigtypeID[itr]);
     s.SetName(trigname);
     s.SetSlotName("trig_pulse_le",1);
     s.SetSignal(trigtime[itr],1);
     s.SetSlotName("trig_pulse_tot",2);
     s.SetSignal(triglength[itr],2);
     trig.AddHit(s);
    } // End of loop over triggers
 
    TMath::Sort(nrtwrevents,twreventtimes,twreventorder,0);
 
    // Add main trigger unless it is already present
    // TODO: select most appropriate trigger as artificial "main"
    // For now: select first trigger in event
    if(!trig.GetHit("main")){
     s.Reset(1);
     s.SetName("main");
     s.SetSlotName("trig_pulse_le",1);
     s.SetSignal(trigtime[0],1);
     s.SetSlotName("trig_pulse_tot",2);
     s.SetSignal(triglength[0],2);
     trig.AddHit(s);
    }
    evt->AddDevice(trig);
 
    //Going to store the MCInfo in the device in case there is any
    if(nmcinfo>0){
      mcinfo.Reset(1);
      for(Int_t imcinfo=0; imcinfo<nmcinfo; imcinfo++){
    TString mcweightitemnamestring(mcweightitemname[imcinfo]);
    mcinfo.AddNamedSlot(mcweightitemname[imcinfo]);
    mcinfo.SetSignal(mcweightitemvalue[imcinfo], mcweightitemname[imcinfo]);
    if(mcweightitemnamestring.Index("OneWeight")>=0) evt->SetWeight(mcweightitemvalue[imcinfo]);
      }
      evt->AddDevice(mcinfo);
    }
    
    // Remember event nr
    lastevent=eventid;
   }
 
   if(waveformtype==1){
    // Amanda module
    omid=aom.GetOMId(String,OM);
    if(omid==681) continue;  // Skip OM 681, which should never give data: avoid all risk of confusion
    om=(IceGOM*)evt->GetIdDevice(omid);
    if (!om)
    {
     aom.Reset(1);
     aom.SetUniqueID(omid);
     evt->AddDevice(aom);
     om=(IceGOM*)evt->GetIdDevice(omid);
    }
    if (!om) continue;
   }
 
   else if(waveformtype==2) {
    // Non-merged waveforms not supported
    cout << "IceRoot::Exec: Non-merged ATWD waveforms not supported in this module." << endl;
    continue;
   }
 
   else if(waveformtype==3 || waveformtype==4) {
    // InIce module
    omid=idom.GetOMId(String,OM);
    om=(IceGOM*)evt->GetIdDevice(omid);
    if (!om)
    {
     idom.Reset(1);
     idom.SetUniqueID(omid);
     evt->AddDevice(idom);
     om=(IceGOM*)evt->GetIdDevice(omid);
    }
    if (!om) continue;
   }
 
   else if(waveformtype==5 || waveformtype==6) {
    // IceTop module
    omid=tdom.GetOMId(String,OM);
    om=(IceGOM*)evt->GetIdDevice(omid);
    if (!om)
    {
     tdom.Reset(1);
     tdom.SetUniqueID(omid);
     evt->AddDevice(tdom);
     om=(IceGOM*)evt->GetIdDevice(omid);
    }
    if (!om) continue;
   }
 
   else {
    cout << "IceRoot::Exec: Unknown waveform type " << waveformtype << endl;
    continue;
   }
 
   // Fill the waveform histogram with this fragment
   // TWR waveform
   if(waveformtype==1){ 
    if(numberbinstwr>0 && twreventnr<nrtwrevents){
     histo.Reset();
     // Store baseline info
     hname="BASELINE-WF";
     hname+=om->GetNwaveforms()+1;
     om->AddNamedSlot(hname);
     om->SetSignal(baseline,hname);
     // Waveform name
     hname="OM";
     hname+=omid;
     hname+="-WF";
     hname+=om->GetNwaveforms()+1;
     hname+="-TWR";
     histo.SetName(hname.Data());
     // Add waveform
     calom=0;
     if(fJEBTDaq) calom=(IceGOM*)fJEBTDaq->GetObject(omid,1);
     if(!calom){
      cout << "IceRoot: No calibration info for OM " << omid << ", skipping this waveform" << endl;
      continue;
     }
     binwidth=calom->GetSignal("BINSIZE");
     if(binwidth<=0){
      cout << "IceRoot: Zero bin width for OM " << omid << ", skipping this waveform" << endl;
      continue;
     }
     twrbuffer=(Int_t)(twrtime/binwidth);
     extstop=(Int_t)calom->GetSignal("EXTSTOP");
     starttime=twreventtimes[twreventorder[twreventnr]]+binwidth*(startbin-twrbuffer+extstop);
     histo.SetBins(numberbinstwr,starttime,starttime+numberbinstwr*binwidth);
     for (Int_t jbin=1; jbin<=numberbinstwr; jbin++)
     {
      histo.SetBinContent(jbin,(Float_t)baseline-(Float_t)twrwaveform[jbin-1]);
     }
     om->SetWaveform(&histo,om->GetNwaveforms()+1);
    } else {
     cout << "IceRoot::Exec: waveformtype=1, but numberbinstwr=0 or nrtwrevents=0." << endl;
    }
   }
 
   // Separate waveforms from DOM launch
   else if(waveformtype==2){
    // Non-merged waveforms not supported
    cout << "IceRoot::Exec: Non-merged ATWD waveform found. You should never get to here anyway." << endl;
    continue;
   }
 
   // Merged ATWD waveform
   else if(waveformtype==3 || waveformtype==5){
    if(numberbinswaveform>0){
     histo.Reset();
     hname="OM";
     hname+=omid;
     hname+="-WF";
     hname+=om->GetNwaveforms()+1;
     hname+="-ATWD";
     histo.SetName(hname.Data());
     starttime=starttimewaveform;
     histo.SetBins(numberbinswaveform,starttime,starttime+numberbinswaveform*binwidth);
     sample.Reset();
     for (Int_t jbin=1; jbin<=numberbinswaveform; jbin++)
     {
      histo.SetBinContent(jbin,waveform[jbin-1]);
     }
     hname="BASELINE-WF";
     hname+=om->GetNwaveforms()+1;
     om->AddNamedSlot(hname);
     om->SetSignal(sample.GetMedian(&histo,2),hname);
     om->SetSignalError(sample.GetSpread(&histo,2),hname);
     om->SetWaveform(&histo,om->GetNwaveforms()+1);
    } else {
     cout << "IceRoot::Exec: waveformtype=" << waveformtype << ", but numberbinswaveform=0." << endl;
    }
   }
 
   // FADC waveform
   else if(waveformtype==4 || waveformtype==6){
    if(numberbinswaveform>0){
     histo.Reset();
     hname="OM";
     hname+=omid;
     hname+="-WF";
     hname+=om->GetNwaveforms()+1;
     hname+="-FADC";
     histo.SetName(hname.Data());
     starttime=starttimewaveform;
     histo.SetBins(numberbinswaveform,starttime,starttime+numberbinswaveform*binwidth);
     for (Int_t jbin=1; jbin<=numberbinswaveform; jbin++)
     {
      histo.SetBinContent(jbin,waveform[jbin-1]);
     }
     om->SetWaveform(&histo,om->GetNwaveforms()+1);
    } else {
     cout << "IceRoot::Exec: waveformtype=" << waveformtype << ", but numberbinswaveform=0." << endl;
    }
   }
 
   // Unknown waveform type
   else {
    cout << "IceRoot::Exec: Unknown waveform type " << waveformtype << endl;
    continue;
   }
 
  } // End of loop over waveforms in input tree
 
  // Write last event to tree
  if(evt->GetDevices("IceGOM"))
  {
   CleanTasks();
   ExecuteTasks(opt);
   if (otree) otree->Fill();
   if (fPrintfreq) { if (!(nevt%fPrintfreq)) evt->HeaderData(); }
   // Update event counter
   nevt++;
  }
 
  // Reset event
  evt->Reset();
 
  // Close input file
  input->Close();
 
  // Clean up
  if (trigtime) delete[] trigtime;
  if (triglength) delete[] triglength;
  if (mcweightitemvalue) delete[] mcweightitemvalue;
  if (twrwaveform) delete[] twrwaveform;
  if (waveform) delete[] waveform;
  if (trackx) delete[] trackx;
  if (tracky) delete[] tracky;
  if (trackz) delete[] trackz;
  if (trackzenith) delete[] trackzenith;
  if (trackazimuth) delete[] trackazimuth;
  if (trackenergy) delete[] trackenergy;
  if (twreventtimes)
  {
   delete[] twreventtimes;
   twreventtimes=0;
  }
  if (twreventorder)
  {
   delete[] twreventorder;
   twreventorder=0;
  }
 
  // Stop looping over input files if max. nr. of events is reached
  if (fMaxevt>-1 && nevt>=fMaxevt)
  {
   break;
  }
 
 } // End of input file loop
 
 // Flush possible memory resident data to the output file
 if (fOutfile) fOutfile->Write();
 
 // Remove the IceEvent object from the environment
 // and delete it as well
 if (evt)
 {
  RemoveObject(evt);
  delete evt;
 }
}