IceF2k.cxx

Go to the documentation of this file.

 
#include "IceF2k.h"
#include "Riostream.h"
 
ClassImp(IceF2k); // Class implementation to enable ROOT I/O
 
IceF2k::IceF2k(const char* name,const char* title) : NcJob(name,title)
{
 
 fSplit=0;
 fBsize=32000;
 fMaxevt=-1;
 fPrintfreq=1;
 fInfiles=0;
 fOutfile=0;
 
 fPdg=0;
 fOmdb=0;
 fFitdefs=0;
 fTrigdefs=0;
 fToffset=0;
 fMctoffset=0;
 fMctracks=3;
 fCompTWR=0;
}

IceF2k::~IceF2k()
{
 
 if (fInfiles)
 {
  delete fInfiles;
  fInfiles=0;
 }
 
 if (fPdg)
 {
  delete fPdg;
  fPdg=0;
 }
 
 if (fOmdb)
 {
  delete fOmdb;
  fOmdb=0;
 }
 
 if (fFitdefs)
 {
  delete fFitdefs;
  fFitdefs=0;
 }
 
 if (fTrigdefs)
 {
  delete fTrigdefs;
  fTrigdefs=0;
 }
}

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

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

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

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

void IceF2k::SetMcToffset(Float_t toffset)
{
 
 fMctoffset=toffset;
}

void IceF2k::SelectMcTracks(Int_t mode)
{
 
 if (mode<0 || mode >3) return;
 fMctracks=mode;
}

void IceF2k::SetCompressedTWR(Int_t mode)
{
 
 if (mode<0 || mode >1) return;
 fCompTWR=mode;
}

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

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

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

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

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

TDatabasePDG* IceF2k::GetPDG()
{
 
 return fPdg;
}

NcObjMatrix* IceF2k::GetOMdbase()
{
 
 return fOmdb;
}

NcDevice* IceF2k::GetFitdefs()
{
 
 return fFitdefs;
}

NcDevice* IceF2k::GetTrigdefs()
{
 
 return fTrigdefs;
}

void IceF2k::Exec(Option_t* opt)
{
 
 if (!fInfiles)
 {
  cout << " *IceF2k Exec* No data input file(s) specified." << endl;
  return;
 }
 
 Int_t ninfiles=fInfiles->GetEntries();
 if (!ninfiles)
 {
  cout << " *IceF2k Exec* No data input file(s) specified." << endl;
  return;
 }
 
 TTree* otree=0;
 if (fOutfile)
 {
  otree=new TTree("T","F2K Data converted to IceEvent structures");
  otree->SetDirectory(fOutfile);
 }
 
 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); 
 
 // Create the particle database and extend it with some F2000 specific definitions
 if (!fPdg) fPdg=new TDatabasePDG();
 Double_t me=fPdg->GetParticle(11)->Mass();
 fPdg->AddParticle("brems"   ,"brems"   ,0,1,0,0,"none",10001001,0,0);
 fPdg->AddParticle("deltae"  ,"deltae"  ,me,1,0,-3,"Lepton",10001002,0,0);
 fPdg->AddParticle("pairprod","pairprod",0,1,0,0,"none",10001003,0,0);
 fPdg->AddParticle("nucl_int","nucl_Int",0,1,0,0,"none",10001004,0,0);
 fPdg->AddParticle("mu_pair" ,"mu_pair" ,0,1,0,0,"none",10001005,0,0);
 fPdg->AddParticle("hadrons" ,"hadrons" ,0,1,0,0,"none",10001006,0,0);
 fPdg->AddParticle("fiberlaser","fiberlaser",0,1,0,0,"none",10002100,0,0);
 fPdg->AddParticle("n2laser"   ,"n2laser"   ,0,1,0,0,"none",10002101,0,0);
 fPdg->AddParticle("yaglaser"  ,"yaglaser"  ,0,1,0,0,"none",10002201,0,0);
 fPdg->AddParticle("z_primary","z_primary",0,1,0,0,"none",10003000,0,0);
 fPdg->AddParticle("a_primary","a_primary",0,1,0,0,"none",10003500,0,0);
 
 // Storage of the used parameters in the IceF2k device
 NcDevice params;
 params.SetNameTitle("IceF2k","IceF2k processor parameters");
 params.SetSlotName("Toffset",1);
 params.SetSlotName("Mctoffset",2);
 params.SetSlotName("Mctracks",3);
 
 // Initialise the job working environment
 SetMainObject(evt);
 if (fOutfile)
 {
  AddObject(fOutfile);
  AddObject(otree);
 }
 
 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 << " F2K 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;
 }
 cout << " TWR input compression level : " << fCompTWR << endl;
 
 ListEnvironment();
 
 // Set DAQ device info
 NcDevice daq;
 daq.SetName("Daq");
 daq.SetSlotName("Muon",1);
 daq.SetSignal(1,1);
 
 Int_t nevt=0;
 Int_t evtsel=0;
 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 input file in the default ascii format (autodetection) for reading 
  fInput=rdmc_mcopen(inputfile.Data(),"r",RDMC_DEFAULT_ASCII_F);
 
  if (!fInput)
  {
   cout << " *IceF2k Exec* No input file found with name : " << inputfile.Data() << endl;
   continue;
  }
 
  // Initialise the event structure 
  rdmc_init_mevt(&fEvent);
 
  // Read the file header information
  rdmc_rarr(fInput,&fHeader);
 
  // Fill the database with geometry, calib. etc... parameters
  // for all the devices
  FillOMdbase();
 
  // Set the fit definitions according to the F2000 header info
  SetFitdefs();
 
  // Set the trigger definitions according to the F2000 header info
  SetTrigdefs();
 
  while (!rdmc_revt(fInput,&fHeader,&fEvent))
  {
   if (fMaxevt>-1 && nevt>=fMaxevt) break;
 
   // Reset the complete Event structure
   evt->Reset();
 
   evt->SetRunNumber(fEvent.nrun);
   evt->SetEventNumber(fEvent.enr);
   evt->SetMJD(fEvent.mjd,fEvent.secs,fEvent.nsecs);
 
   evt->AddDevice(daq);
 
   // Take trigger offset into account which might have been
   // introduced during the filtering process.
   // For simulated data this will be treated separately in PutMcTracks().
   fToffset=fEvent.t_offset;
 
   PutTrigger();
 
   PutMcTracks();
 
   PutRecoTracks();
 
   PutHits();
 
   // Enter the IceF2k processor parameters into the event structure 
   params.SetSignal(fToffset,1);
   params.SetSignal(fMctoffset,2);
   params.SetSignal(fMctracks,3);
   evt->AddDevice(params);
 
   // Invoke all available sub-tasks (if any)
   CleanTasks();
   ExecuteTasks(opt);
 
   if (fPrintfreq)
   {
    if (!(nevt%fPrintfreq)) evt->HeaderData();
   }
 
   // Write the complete structure to the output Tree for accepted events
   evtsel=1;
   NcDevice* seldev=(NcDevice*)evt->GetDevice("NcEventSelector");
   if (seldev)
   {
    if (seldev->GetSignal("Select") < 0.1) evtsel=-1;
   }
   if (otree && evtsel==1) otree->Fill();
 
   // Update event counter
   nevt++;
  }
  if (fMaxevt>-1 && nevt>=fMaxevt) break;
 }
 
 // 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;
 }
}

void IceF2k::FillOMdbase()
{
 
 if (fHeader.nch<=0)
 {
  if (fOmdb)
  {
   delete fOmdb;
   fOmdb=0;
  }
  return;
 }
 
 Int_t adccal=fHeader.is_calib.adc;
 Int_t tdccal=fHeader.is_calib.tdc;
 Int_t totcal=fHeader.is_calib.tot;
 
 TF1 fadccal("fadccal","(x-[1])*[0]");
 TF1 fadcdecal("fadcdecal","(x/[0])+[1]");
 fadccal.SetParName(0,"BETA-ADC");
 fadccal.SetParName(1,"PED-ADC");
 fadcdecal.SetParName(0,"BETA-ADC");
 fadcdecal.SetParName(1,"PED-ADC");
 
 TF1 ftdccal("ftdccal","(x*[0])-[1]-([0]-1.)*32767.-[2]/sqrt([3])");
 TF1 ftdcdecal("ftdcdecal","(x+([0]-1.)*32767.+[1]+[2]/sqrt([3]))/[0]");
 ftdccal.SetParName(0,"BETA-TDC");
 ftdccal.SetParName(1,"T0");
 ftdccal.SetParName(2,"ALPHA-TDC");
 ftdccal.SetParName(3,"ADC-SLEW");
 ftdcdecal.SetParName(0,"BETA-TDC");
 ftdcdecal.SetParName(1,"T0");
 ftdcdecal.SetParName(2,"ALPHA-TDC");
 ftdcdecal.SetParName(3,"ADC-SLEW");
 
 TF1 ftotcal("ftotcal","x*[0]");
 TF1 ftotdecal("ftotdecal","x/[0]");
 ftotcal.SetParName(0,"BETA-TOT");
 ftotdecal.SetParName(0,"BETA-TOT");
 
 if (fOmdb)
 {
  fOmdb->Reset();
 }
 else
 {
  fOmdb=new NcObjMatrix();
  fOmdb->SetNameTitle("OMDBASE","The OM geometry, calib. etc... database");
  fOmdb->SetOwner();
 }
 
 IceAOM* dev=0;
 Double_t pos[3]={0,0,0};
 for (Int_t i=0; i<fHeader.nch; i++)
 {
  dev=new IceAOM();
  dev->SetUniqueID(i+1);
  // Slots to hold the various (de)calibration functions  
  dev->SetSlotName("ADC",1);
  dev->SetSlotName("LE",2);
  dev->SetSlotName("TOT",3);
  // Slots to hold hardware parameters
  dev->SetSlotName("TYPE",4);
  dev->SetSlotName("ORIENT",5);
  dev->SetSlotName("THRESH",6);
  dev->SetSlotName("SENSIT",7);
  dev->SetSlotName("READOUT",8); // 0=unknown 1=electrical 2=optical 3=digital
 
  pos[0]=fHeader.x[i];
  pos[1]=fHeader.y[i];
  pos[2]=fHeader.z[i];
  dev->SetPosition(pos,"car");
 
  fadccal.SetParameter(0,fHeader.cal[i].beta_a);
  fadccal.SetParameter(1,fHeader.cal[i].ped);
  fadcdecal.SetParameter(0,fHeader.cal[i].beta_a);
  if (!fHeader.cal[i].beta_a) fadcdecal.SetParameter(0,1);
  fadcdecal.SetParameter(1,fHeader.cal[i].ped);
 
  ftdccal.SetParameter(0,fHeader.cal[i].beta_t);
  ftdccal.SetParameter(1,fHeader.cal[i].t_0);
  ftdccal.SetParameter(2,fHeader.cal[i].alpha_t);
  ftdccal.SetParameter(3,1.e20);
  ftdcdecal.SetParameter(0,fHeader.cal[i].beta_t);
  if (!fHeader.cal[i].beta_t) ftdcdecal.SetParameter(0,1);
  ftdcdecal.SetParameter(1,fHeader.cal[i].t_0);
  ftdcdecal.SetParameter(2,fHeader.cal[i].alpha_t);
  ftdcdecal.SetParameter(3,1.e20);
 
  ftotcal.SetParameter(0,fHeader.cal[i].beta_tot);
  ftotdecal.SetParameter(0,fHeader.cal[i].beta_tot);
  if (!fHeader.cal[i].beta_tot) ftotdecal.SetParameter(0,1);
 
  if (adccal)
  {
   dev->SetDecalFunction(&fadcdecal,1);
  }
  else
  {
   dev->SetCalFunction(&fadccal,1);
  }
 
  if (tdccal)
  {
   dev->SetDecalFunction(&ftdcdecal,2);
  }
  else
  {
   dev->SetCalFunction(&ftdccal,2);
  }
 
  if (totcal)
  {
   dev->SetDecalFunction(&ftotdecal,3);
  }
  else
  {
   dev->SetCalFunction(&ftotcal,3);
  }
 
  dev->SetSignal(fHeader.type[i],4);
  dev->SetSignal((Float_t)fHeader.costh[i],5);
  dev->SetSignal(fHeader.thresh[i],6);
  dev->SetSignal(fHeader.sensit[i],7);
  dev->SetSignal(0.,8);
 
  fOmdb->EnterObject(i+1,1,dev);
 }
}

void IceF2k::SetFitdefs()
{
 
 if (fHeader.n_fit<=0)
 {
  if (fFitdefs)
  {
   delete fFitdefs;
   fFitdefs=0;
  }
  return;
 }
 
 if (fFitdefs)
 {
  fFitdefs->Reset(1);
 }
 else
 {
  fFitdefs=new NcDevice();
 }
 
 fFitdefs->SetName("FitDefinitions");
 fFitdefs->SetHitCopy (1);
 
 NcSignal s;
 s.Reset();
 
 for (Int_t i=0; i<fHeader.n_fit; i++)
 {
  s.SetUniqueID(fHeader.def_fit[i].id);
  s.SetName(TString(fHeader.def_fit[i].tag));
 
  for (Int_t j=0; j<fHeader.def_fit[i].nwords; j++)
  {
   s.SetSlotName(TString(fHeader.def_fit[i].words[j]),j+1);
   s.SetSignal(0,j+1);
  }
 
  fFitdefs->AddHit(s);
  s.Reset(1);
 }
}

void IceF2k::SetTrigdefs()
{
 
 if (fHeader.n_trigger<=0)
 {
  if (fTrigdefs)
  {
   delete fTrigdefs;
   fTrigdefs=0;
  }
  return;
 }
 
 if (fTrigdefs)
 {
  fTrigdefs->Reset(1);
 }
 else
 {
  fTrigdefs=new NcDevice();
 }
 
 fTrigdefs->SetName("TrigDefinitions");
 fTrigdefs->SetHitCopy (1);
 
 NcSignal s;
 s.Reset();
 
 for (Int_t i=0; i<fHeader.n_trigger; i++)
 {
  s.SetUniqueID(fHeader.def_trig[i].id);
  s.SetName(TString(fHeader.def_trig[i].tag));
 
  for (Int_t j=0; j<fHeader.def_trig[i].nwords; j++)
  {
   s.SetSlotName(TString(fHeader.def_trig[i].words[j]),j+1);
   s.SetSignal(0,j+1);
  }
 
  fTrigdefs->AddHit(s);
  s.Reset(1);
 }
}

void IceF2k::PutMcTracks()
{
 
 IceEvent* evt=(IceEvent*)GetMainObject();
 if (!evt || fEvent.ntrack<=0) return;
 
 // User defined trigger offset in case of simulated data.
 // The offset in the F2K file is meant to put the primary interaction
 // well ahead of the detector trigger.
 // See the introductory docs of this IceF2k class for further details.
 fToffset=fMctoffset;
 
 if (!fMctracks) return;
 
 // Loop over all the tracks and add them to the current event
 NcTrack t;
 Double_t vec[3];
 NcPosition r;
 Nc3Vector p;
 Int_t tid=0;
 Int_t idpdg=0;
 Int_t idf2k=0;
 for (Int_t i=0; i<fEvent.ntrack; i++)
 {
  t.Reset ();
 
  // Beginpoint of the track
  vec[0]=fEvent.gen[i].x;
  vec[1]=fEvent.gen[i].y;
  vec[2]=fEvent.gen[i].z;
  r.SetPosition(vec,"car");
  t.SetBeginPoint(r);
 
  // Endpoint of the track
  vec[0]+=fEvent.gen[i].length*fEvent.gen[i].px;
  vec[1]+=fEvent.gen[i].length*fEvent.gen[i].py;
  vec[2]+=fEvent.gen[i].length*fEvent.gen[i].pz;
  r.SetPosition(vec,"car");
  t.SetEndPoint(r);
 
  // Momentum in GeV/c
  vec[0]=fEvent.gen[i].e*fEvent.gen[i].px*1e-3;
  vec[1]=fEvent.gen[i].e*fEvent.gen[i].py*1e-3;
  vec[2]=fEvent.gen[i].e*fEvent.gen[i].pz*1e-3;
  p.SetVector (vec,"car");
  t.Set3Momentum(p);
 
  // MC tracks are indicated by negative track id's
  tid=fEvent.gen[i].tag;
  t.SetId(-abs(tid));
 
  idf2k=fEvent.gen[i].id;
  idpdg=0;
  if (idf2k>1000)
  {
   idpdg=idf2k+10000000;
  }
  else if (idf2k <= 48)
  {
   idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
  }
  else
  {
   if (idf2k==201) idpdg=12;
   if (idf2k==202) idpdg=14;
   if (idf2k==203) idpdg=16;
   if (idf2k==204) idpdg=-12;
   if (idf2k==205) idpdg=-14;
   if (idf2k==206) idpdg=-16;
  }
 
  // Check for the user selected MC track storage
  if (fMctracks==1) // Store only muon and muon-neutrino tracks
  {
   if (abs(idpdg)!=13 && abs(idpdg)!=14) continue;
  }
  else if (fMctracks==2) // Store all lepton tracks
  {
   if (abs(idpdg)<11 || abs(idpdg)>16) continue;
  }
 
  t.SetParticleCode(idpdg);
  t.SetName(fPdg->GetParticle(idpdg)->GetName());
  t.SetTitle("MC track");
  t.SetMass(fPdg->GetParticle(idpdg)->Mass());
  t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
 
  evt->AddTrack(t);
 }
 
 // Create the pointers to each particle's parent particle.
 Int_t txid=0;
 Int_t parid=0;
 for (Int_t itk=1; itk<=evt->GetNtracks (); itk++)
 {
  NcTrack* tx=evt->GetTrack(itk);
 
  if (!tx) continue;
 
  txid=tx->GetId();
 
  parid=-1;
  for (Int_t j=0; j<fEvent.ntrack; j++)
  {
   tid=fEvent.gen[j].tag;
   if (-abs(tid) == txid) parid=fEvent.gen[j].parent;
  }
 
  if (parid<0) continue;
 
  NcTrack* tpar=evt->GetIdTrack(-abs(parid));
 
  if (!tpar) continue;
 
  tx->SetParentTrack(tpar);
 }
}

void IceF2k::PutRecoTracks()
{
 
 IceEvent* evt=(IceEvent*)GetMainObject();
 if (!evt || fEvent.nfit<=0) return;
 
 // Loop over all the tracks and add them to the current event
 NcTrack t;
 Double_t vec[3];
 NcPosition r;
 Nc3Vector p;
 Int_t tid=0;
 Int_t idpdg=0;
 Int_t idf2k=0;
 for (Int_t i=0; i<fEvent.nfit; i++)
 {
  t.Reset ();
 
  // Beginpoint of the track
  vec[0]=fEvent.rec[i].x;
  vec[1]=fEvent.rec[i].y;
  vec[2]=fEvent.rec[i].z;
  r.SetPosition(vec,"car");
  t.SetBeginPoint(r);
 
  // Endpoint of the track
  vec[0]+=fEvent.rec[i].length*fEvent.rec[i].px;
  vec[1]+=fEvent.rec[i].length*fEvent.rec[i].py;
  vec[2]+=fEvent.rec[i].length*fEvent.rec[i].pz;
  r.SetPosition(vec,"car");
  t.SetEndPoint(r);
 
  // Momentum in GeV/c
  if (fEvent.rec[i].e > 0)
  {
   vec[0]=fEvent.rec[i].e*fEvent.rec[i].px*1e-3;
   vec[1]=fEvent.rec[i].e*fEvent.rec[i].py*1e-3;
   vec[2]=fEvent.rec[i].e*fEvent.rec[i].pz*1e-3;
  }
  else // Give the track a nominal momentum of 1 GeV/c
  {
   vec[0]=fEvent.rec[i].px;
   vec[1]=fEvent.rec[i].py;
   vec[2]=fEvent.rec[i].pz;
  }
  p.SetVector (vec,"car");
  t.Set3Momentum(p);
 
  // Use the fit number as track id
  tid=fEvent.rec[i].tag;
  t.SetId(abs(tid));
 
  idf2k=fEvent.rec[i].id;
  idpdg=0;
  if (idf2k>1000)
  {
   idpdg=idf2k+10000000;
  }
  else if (idf2k <= 48)
  {
   idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
  }
  else
  {
   if (idf2k==201) idpdg=12;
   if (idf2k==202) idpdg=14;
   if (idf2k==203) idpdg=16;
   if (idf2k==204) idpdg=-12;
   if (idf2k==205) idpdg=-14;
   if (idf2k==206) idpdg=-16;
  }
 
  t.SetParticleCode(idpdg);
  t.SetNameTitle("Sieglinde","RECO track");
  t.SetMass(fPdg->GetParticle(idpdg)->Mass());
  t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
 
  // Retrieve the various fit parameters for this track
  NcSignal* fitdata=fFitdefs->GetIdHit(i);
  for (Int_t jval=0; jval<fEvent.fresult[i].nval; jval++)
  {
   fitdata->SetSignal(fEvent.fresult[i].val[jval],jval+1);
  }
 
  // Store the various fit parameters for this track
  t.SetFitDetails(fitdata);
 
  // Store the various reco tracks as track hypotheses.
  // A copy of the first reco track is entered as a new track instance
  // into the event and all reco tracks (incl. the first one) are
  // stored as hypotheses linked to this new reco track.
  if (i==0)
  {
   evt->AddTrack(t);
   NcTrack* tx=evt->GetTrack(evt->GetNtracks());
   Int_t nrec=evt->GetNtracks(1);
   tx->SetId(nrec+1);
  }
  NcTrack* tx=evt->GetTrack(evt->GetNtracks());
  if (tx) tx->AddTrackHypothesis(t);
 }
}

void IceF2k::PutHits()
{
 
 IceEvent* evt=(IceEvent*)GetMainObject();
 if (!evt) return;
 
 // Loop over all the hits and add them to the current event
 IceAOM om;
 NcSignal s;
 s.SetSlotName("ADC",1);
 s.SetSlotName("LE",2);
 s.SetSlotName("TOT",3);
 Int_t chan=0;
 Int_t maxchan=800;
 if (fOmdb) maxchan=fHeader.nch;
 IceAOM* omx=0;
 NcSignal* sx=0;
 Int_t tid=0;
 NcTrack* tx=0;
 Float_t adc=0;
 Float_t adcfirst=0; // Adc value of the first hit of an OM
 for (Int_t i=0; i<fEvent.nhits; i++)
 {
  chan=fEvent.h[i].ch+1;
  if (chan>maxchan) continue; // Channels 9001, 9002 etc are trigger channels
 
  // Get corresponding device from the current event structure  
  omx=(IceAOM*)evt->GetIdDevice(chan);
  if (!omx)
  {
   if (fOmdb)
   {
    omx=(IceAOM*)fOmdb->GetObject(chan,1);
    evt->AddDevice(omx);
   }
   else
   {
    om.Reset(1);
    om.SetUniqueID(chan);
    evt->AddDevice(om);
   }
   omx=(IceAOM*)evt->GetIdDevice(chan);
  }
 
  if (!omx) continue;
 
  adc=fEvent.h[i].amp;
 
  // Multiple hits in the same OM with the same ADC value
  // are indicated by "*" in the F2K file.
  // This corresponds to a value of -2 in the data structure.
  if (int(adc) == -2)
  {
   adc=adcfirst;
  }
  else
  {
   adcfirst=adc;
  }
  s.Reset();
  s.SetUniqueID(fEvent.h[i].id);
  s.SetSignal(adc,1);
  s.SetSignal((fEvent.h[i].t-fToffset),2);
  s.SetSignal(fEvent.h[i].tot,3);
 
  omx->AddHit(s);
 
  sx=omx->GetHit(omx->GetNhits());
  if (!sx) continue;
 
  // ADC dependent TDC (de)calibration function for this hit
  TF1* fcal=omx->GetCalFunction("LE");
  TF1* fdecal=omx->GetDecalFunction("LE");
  if (fcal) sx->SetCalFunction(fcal,2);
  if (fdecal) sx->SetDecalFunction(fdecal,2);
  fcal=sx->GetCalFunction(2);
  fdecal=sx->GetDecalFunction(2);
  adc=sx->GetSignal(1,-4);
  if (adc>0)
  {
   if (fcal) fcal->SetParameter(3,adc);
   if (fdecal) fdecal->SetParameter(3,adc);
  }
  else
  {
   if (fcal) fcal->SetParameter(3,1.e20);
   if (fdecal) fdecal->SetParameter(3,1.e20);
  }
 
  // Bi-directional link between this hit and the track that caused the ADC value.
  // This F2K info is probably only present for MC tracks.
  tid=fEvent.h[i].ma;
  if (tid > 0)
  {
   tx=evt->GetIdTrack(tid); // Reco tracks
   if (!tx) tx=evt->GetIdTrack(-tid); // MC tracks
   if (tx) sx->AddTrack(*tx);
  }
  else
  {
   if (tid == -2) sx->SetNameTitle("N","Noise");
   if (tid == -3) sx->SetNameTitle("A","Afterpulse");
  }
 }
 
 // Loop over all the waveforms and add the histo(s) to the corresponding OM's
 TH1F histo;
 Int_t nbins=0;
 Float_t xlow=0;
 Float_t xup=0;
 TString hname;
 Float_t amp=0; // Amplitude sampler for compressed waveforms
 for (Int_t iwf=0; iwf<fEvent.nwf; iwf++)
 {
  chan=fEvent.wf[iwf].om;
  if (chan<=0 || chan>maxchan) continue; // Skip trigger channels
 
  // Get corresponding device from the current event structure  
  omx=(IceAOM*)evt->GetIdDevice(chan);
  if (!omx)
  {
   if (fOmdb)
   {
    omx=(IceAOM*)fOmdb->GetObject(chan,1);
    evt->AddDevice(omx);
   }
   else
   {
    om.Reset(1);
    om.SetUniqueID(chan);
    evt->AddDevice(om);
   }
   omx=(IceAOM*)evt->GetIdDevice(chan);
  }
 
  if (!omx) continue;
 
  hname="BASELINE-WF";
  hname+=omx->GetNwaveforms()+1;
  omx->AddNamedSlot(hname);
  omx->SetSignal(fEvent.wf[iwf].baseline,hname);
 
  // Fill the waveform histogram
  hname="OM";
  hname+=chan;
  hname+="-WF";
  hname+=omx->GetNwaveforms()+1;
 
  histo.Reset();
  histo.SetName(hname.Data());
  nbins=fEvent.wf[iwf].ndigi;
  xlow=fEvent.wf[iwf].t_start;
  xup=xlow+float(nbins)*fEvent.wf[iwf].t_bin;
  histo.SetBins(nbins,xlow,xup);
 
  for (Int_t jbin=1; jbin<=fEvent.wf[iwf].ndigi; jbin++)
  {
   if (!fCompTWR) // Uncompressed waveform input
   {
    histo.SetBinContent(jbin,fEvent.wf[iwf].baseline-fEvent.wf[iwf].digi[jbin-1]);
   }
   else // Compressed waveform input
   {
    if (jbin==1) amp=fEvent.wf[iwf].digi[0];
    if (jbin==2) amp=fEvent.wf[iwf].digi[1]-histo.GetBinContent(1);
    if (jbin>2) amp=fEvent.wf[iwf].digi[jbin-1]-2.*histo.GetBinContent(jbin-1)+histo.GetBinContent(jbin-2);
    histo.SetBinContent(jbin,-amp);
   }
  }
 
  omx->SetWaveform(&histo,omx->GetNwaveforms()+1);
 }
 
 // Set bi-directional links between hits and reco track hypotheses.
 // Note : Reco tracks are recognised by their positive id.
 Int_t hid=0;
 TObjArray* rectracks=evt->GetTracks(1);
 for (Int_t jtk=0; jtk<rectracks->GetEntries(); jtk++)
 {
  tx=(NcTrack*)rectracks->At(jtk);
  if (!tx) continue;
  
  for (Int_t jhyp=1; jhyp<=tx->GetNhypotheses(); jhyp++)
  {
   NcTrack* hypx=tx->GetTrackHypothesis(jhyp);
   if (!hypx) continue;
 
   // Loop over all combinations of F2K fits and used OM hits
   for (Int_t k=0; k<fEvent.nfit_uses; k++)
   {
    if (fEvent.fit_uses[k].useid != hypx->GetId()) continue;
    hid=fEvent.fit_uses[k].hitid;
    sx=evt->GetIdHit(hid,"IceAOM");
    if (sx) sx->AddTrack(*hypx);
   }
  }
 }
}

void IceF2k::PutTrigger()
{
 
 if (!fTrigdefs) return;
 
 IceEvent* evt=(IceEvent*)GetMainObject();
 if (!evt || fEvent.ntrig<=0) return;
 
 NcDevice trig;
 trig.SetNameTitle("Trigger","Amanda/IceCube event triggers");
 NcSignal s;
 TString trigname;
 TString slotname;
 Int_t id=0;
 Int_t nval=0;
 for (Int_t i=0; i<fEvent.ntrig; i++)
 {
  id=fEvent.ptrig[i].id;
  nval=fEvent.ptrig[i].nval;
  if (!nval) continue;
  NcSignal* tdef=fTrigdefs->GetIdHit(id+1);
  if (!tdef) continue;
  trigname=tdef->GetName();
  s.Reset(1);
  s.SetName(trigname);
  s.SetUniqueID(id+1);
  for (Int_t jval=0; jval<fEvent.ptrig[i].nval; jval++)
  {
   slotname=tdef->GetSlotName(jval+1);
   s.SetSlotName(slotname,jval+1);
   s.SetSignal(fEvent.ptrig[i].val[jval],jval+1);
  }
  trig.AddHit(s);
 }
 
 // Store the trigger data into the IceEvent structure
 evt->AddDevice(trig);
}